В заметке Тепловое взаимодействие ячеек в соте было промоделировано изменение температуры для двустороннего кластера сот размером 11 на 11 ячеек при помещении туда четырех пчел и исходной температуре 12°С. Для каждой пустой ячейки решалась система из 10 уравнений - по количеству граней плюс уравнение для воздуха, а для ячейки с пчелой - еще дополнительно 3 уравнения: для торакса, головы и брюшка. Решение находилось методом fsolve библиотеки scipy - это численный метод (обычно алгоритм Минпака), который внутри себя подбирает корни системы нелинейных уравнений. Моделирование производилось с небольшим шагом по времени - 0.01 с, т.к. увеличение шага приводило к расхождению решения (примерно с шага 0.03 с). Поскольку для достижения требуемой точности на каждом шаге проводилось множество итераций, решение находилось медленно. Десять минут (шестьсот секунд) реального времени моделировалось примерно за 300 секунд. Это значит, что рамка Дадана считалась бы около часа, а весь улей более десяти часов. Десять минут - за 10 часов! Это, конечно же, неприемлемо, если иметь в виду цель - создание цифрового двойника семьи в улье.
Я попытался где можно векторизовать уравнения и применить известные методы ускорения в Python, в частности надстройку numba, но это ускорило вычисления не более, чем в полтора раза. Тогда я опять обратился за помощью к ИИ, и на этот раз он справился (мой первый опыт, как вы помните, был не очень удачным).
Суть предложений ИИ сводилась к двум: изменению алгоритма решения задачи (системы уравнений) на метод Эйлера и жесткой структуре данных в виде массивов: массива температур, массива топологии (связей), массива теплоемкостей и т.д. Оба эти предложения нацелены на то, чтобы дать надстройке numba проявить свою максимальную эффективность
По алгоритму
Сопоставим теплоемкости:
- Теплоемкость воска ячейки \(C_{wax}\) : 124 мДж/°C.
- Теплоемкость пчелы \(C_{bee}\): 349 мДж/°C.
- Теплоемкость воздуха в ячейке с пчелой \(C_{air\_bee}\): 0.145 мДж/°C.
- Теплоемкость воздуха в пустой ячейке \(C_{air\_cell}\): около 0.367 мДж/°C.
Разница между воздухом и воском — в сотни раз. Воздух внутри микро-объема ячейки является малоинерционной средой. Поэтому, можно считать воздух в ячейке безынерционным - пренебречь его теплоемкостью. В этом случае температура воздуха в ячейке на текущем шаге — это просто средневзвешенная температура окружения:
$$T_{ai}=\frac{\frac{T_{th}}{R_{th\_ai}}+\frac{T_{hd}}{R_{hd\_ai}}+\frac{T_{ab}}{R_{ab\_ai}}+\sum _{k=0}^{8}\frac{T_{wall}[k]}{R_{ai\_wall}[k]}+\frac{T_{ambient}}{R_{ab\_}}}{\frac{1}{R_{th\_ai}}+\frac{1}{R_{hd\_ai}}+\frac{1}{R_{ab\_ai}}+\sum _{k=0}^{8}\frac{1}{R_{ai\_wall}[k]}+\frac{1}{R_{ab\_}}}$$
Прямой пересчет температур (явная схема Эйлера).
Раз у нас есть \(T_{ai}\), новые температуры для следующего шага времени (t + dt) вычисляются за 1 такт безо всяких итераций, например, для температуры головы:
$$T_{hd}=T_{hd\_prev}+\frac{dt}{C_{hd}}\cdot \left(p_{s\_hd\_prev}+\frac{T_{th\_prev}-T_{hd\_prev}}{R_{th\_hd}}-\frac{T_{hd\_prev}-T_{ai\_prev}}{R_{hd\_ai}}-q_{evap\_hd\_prev}\right)$$
Аналогично для \(T_{hd}, T_{th}\) и всех 9 граней \(T_{wall}\).
По массивам
Основные массивы и их формы:
- Массив теплоемкостей
C_grid (N_cells×13) - Массив температур ячеек
T_grid (N_cells×13) - Массив топологии
Topology (N_cells×9) - Массив соседей по граням
Neighbor_Faces(N_cells×9) - Массив статусов
States(N_cells)
Такое структурирование данных позволяет при моделировании одного шага во времени убрать цикл по ячейкам, не перебирать составляющие (13 - девять стенок, три части тела пчелы и воздух). Т.е. в программе остается один цикл - по времени.
Эти предложения сам же ИИ и оформил в скрипт.
Задав такой же кластер 11×11 ячеек и поместив там 4 пчелы, можно произвести сравнение результатов нового алгоритма и старого. Напомню, ранее такой расчет для 600 секунд реального времени выполнялся примерно за 300 секунд при шаге моделирования по времени в 0.01 с, т.е. число шагов составляло 60000,а время выполнения одного шага - 0.05 с. При таком же шаге моделирования новая программа выполнила расчет этих 600 с примерно за 0.8 секунды, т.е. в 375 раз быстрее!
Этот результат уже может позволить моделировать улей быстрее реального времени, и вновь появляется перспектива создания цифрового двойника. Однако, что насчет устойчивости? Мы отмечали, что при увеличении шага примерно с 0.03 секунд устойчивость теряется.
Постепенное увеличение шага в новом алгоритме показало, что устойчивость теряется примерно с шага в 0.3 секунды. Это значит, что можно с запасом поставить шаг моделирования 0.1 секунды, т.е. увеличить его в 10 раз по сравнению со старым алгоритмом. Общее время моделирования 600 секунд для нашего кластера 11×11 ячеек в этом случае составило 0.08 секунды, что в 3750 раз быстрее исходного.
В этот момент чат с ИИ практически переполнился контекстом (это особенность современного этапа развития всех ИИ-чатов), а текст программы в ходе доработки стал не очень читабельным. Поэтому я попросил другой ИИ-чат подготовить текст программы к публикации. "Свежий" ИИ творчески отнесся к делу и предложил еще две возможности поднять быстродействие:
- добавил
fastmath=Trueв обертку@njit(ускоряет математические операции). - убрал
.copy()внутри шага симуляции; вместо этого используется классический трюк "двойного буфера" (double-buffering): мы передаем два массива (T_currentиT_next), пишем вT_next, а затем просто меняем их местами:T_cur,T_next=T_next, T_cur.Это убирает аллокацию памяти на каждом шаге и снимает нагрузку с Garbage Collector.
В результате время моделирования уменьшилось до 0.055 с, т.е. более, чем в 5000 раз по сравнению с исходным.
Быстродействие, таким образом, увеличили. Но известно, что метод Эйлера сам по себе вносит погрешность в решение для нелинейных систем (а она у нас существенно нелинейная). Кроме того, в новом алгоритме мы пренебрегли теплоемкостью воздуха при расчете температуры воздуха в ячейке, что тоже могло внести дополнительную погрешность. Сравним алгоритмы по точности.
 |
 |
 |
 |
 |
 |
В левой колонке результаты старого алгоритма, в правой - нового. Видим, что качественно картина не поменялась - пчелы греют сот вокруг себя. Но расхождения действительно есть: в новом алгоритме мощность, отдаваемая пчелами несколько больше, поэтому и температуры в ячейках с пчелами и окрестностях - выше. При этом разница составляет менее 0.2°С. Это является платой за повышение быстродействия. На данном этапе полагаю такую погрешность допустимой для моделирования. Если в дальнейшем окажется, что такой точности недостаточно - займемся повышением точности, будем решать проблемы по мере их поступления.
P.S. Кстати рамка Дадана, которая ранее считалась бы более часа, сейчас считается менее двух секунд:
🐝 Python-скрипт
""" Программа к заметке "Повышение производительности" """
import time
import math
import numpy as np
from numba import njit
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
# ==========================================================
# 1. ФИЗИЧЕСКИЕ ПАРАМЕТРЫ И ГЕОМЕТРИЯ (Бывший config_11x11.py)
# ==========================================================
# Сетка и условия
cols, rows = 11, 11
variant = 2 # 1: Тонкий сот, 2: Коричневый сот, 3: Толстый сот
T_ambient = 12.0 # Начальная температура воздуха (°C)
# Пчела
d_bee = 0.004
s_r = math.pi * (d_bee / 2) ** 2
h_head, h_thorax, h_ab = 0.0015, 0.0045, 0.006
h_total = h_head + h_thorax + h_ab
A_hd = math.pi * d_bee * h_head + s_r
A_th = math.pi * d_bee * h_thorax
A_ab = math.pi * d_bee * h_ab + s_r
V_bee = h_total * s_r
ke0 = 3e-5
Pmin = 0.2e-3
k_hd, k_ab = 0.065, 0.042
m_head, m_thorax, m_ab = 10.25e-6, 32.5e-6, 57e-6
c = 3500
C_hd, C_th, C_ab = m_head * c, m_thorax * c, m_ab * c
# Воздух
ρ_air, k_air, c_air = 1.205, 0.025, 1005
# Воск и ячейки
L_w = [0.0001, 0.00014, 0.0002]
L_b = [0.00013, 0.00055, 0.0012]
c_w = [3300., 3000., 2700.]
λ_w = [0.235, 0.1608, 0.055]
eps_w = [0.22, 0.5, 0.9]
L_wall = L_w[variant - 1]
L_bottom = L_b[variant - 1]
c_wax = c_w[variant - 1]
λ_wax = λ_w[variant - 1]
ρ_wax = 960
thread = 5.4e-3
d_cell = thread - L_wall
h_cell = 12e-3
dh = 1e-3
a = d_cell / 2. / np.cos(np.radians(30))
A_lateral = a * (h_cell - dh / 2.)
d1, d2 = d_cell, np.sqrt(a ** 2 + dh ** 2)
A_bottom = d1 * d2 / 2.
C_wax_lateral = ρ_wax * A_lateral * L_wall * c_wax
C_wax_bottom = ρ_wax * A_bottom * L_bottom * c_wax
A_wall_total = 6 * A_lateral + 3 * A_bottom
area_hex = (3 * math.sqrt(3) / 2) * a ** 2
A_entrance = area_hex
V_cell = area_hex * h_total
C_air_cell = ρ_air * V_cell * c_air
C_air_bee = ρ_air * (V_cell - V_bee) * c_air
# Излучение и тепловые сопротивления
sigma = 5.67e-8
t0 = -273.15
eps = 0.22
eps_wall = eps_w[variant - 1]
se = sigma * eps_wall
se_th, se_hd, se_ab_ = se * A_th, se * (A_hd - s_r), se * (A_ab - s_r)
h_in = k_air / ((d_cell - d_bee) / 2)
h_entrance = 11.0
R_th_ai = 1 / (h_in * A_th)
R_hd_ai = 1 / (h_in * A_hd)
R_ab_ai = 1 / (h_in * (A_ab - s_r))
R_ab_ = 1 / (h_entrance * s_r)
hA_ = h_entrance * (A_entrance - s_r)
hA = h_entrance * A_entrance
R_head_air_open = 1 / (h_entrance * A_hd)
R_thorax_air_open = 1 / (h_entrance * A_th)
R_ab_air_open = 1 / (h_entrance * A_ab)
# Режимы работы (0.0 - обычный, 45e-3 - режим печки)
P_warm_up = 0.0
# ==========================================================
# 2. ГЕНЕРАЦИЯ ТОПОЛОГИИ И СЕТКИ
# ==========================================================
CELL_EMPTY, CELL_BEE = 0, 5
I_HD, I_TH, I_AB, I_AI = 0, 1, 2, 3
# Массивы сопротивлений и теплоемкостей для Numba
R_ai_wall_arr = np.array([L_wall / (k_air * A) for A in [A_lateral] * 6 + [A_bottom] * 3], dtype=np.float64)
R_wall_self = np.array([L_wall / (λ_wax * A_lateral)] * 6 + [L_bottom / (λ_wax * A_bottom)] * 3, dtype=np.float64)
C_wall_base = np.array([C_wax_lateral] * 6 + [C_wax_bottom] * 3, dtype=np.float64)
# Генерация координат сот
cells = []
for side in ["Y", "^"]:
for col in range(-int(cols / 2), int(cols / 2) + 1):
for row in range(-int(rows / 2), int(rows / 2) + 1):
cells.append((0, side, col, row))
cell_to_index = {cell: i for i, cell in enumerate(cells)}
N_CELLS = len(cells)
EVEN_ROW_DELTAS = [(-1, 0), (-1, -1), (0, -1), (1, 0), (0, 1), (-1, 1)]
ODD_ROW_DELTAS = [(1, 0), (1, -1), (0, -1), (-1, 0), (0, 1), (1, 1)]
def get_neighbor_for_face(side, col, row, face_id):
if face_id < 6:
deltas = ODD_ROW_DELTAS if (row % 2 == 1) else EVEN_ROW_DELTAS
if side == "^":
deltas = EVEN_ROW_DELTAS if (row % 2 == 1) else ODD_ROW_DELTAS
dc, dr = deltas[face_id]
return (0, side, int(col + dc), int(row + dr))
else:
k = face_id - 6
if side == "Y":
bottom_list = [(col, row - 1), (col, row), (col - 1, row)] if row % 2 == 0 else [(col, row - 1), (col, row), (col + 1, row)]
n_side = "^"
else:
bottom_list = [(col, row), (col, row + 1), (col + 1, row)] if row % 2 == 0 else [(col, row), (col, row + 1), (col - 1, row)]
n_side = "Y"
n_col, n_row = bottom_list[k]
return (0, n_side, n_col, n_row)
# Построение матриц топологии
Topology = np.full((N_CELLS, 9), -1, dtype=np.int32)
Neighbor_Faces = np.full((N_CELLS, 9), -1, dtype=np.int32)
all_cells_set = set(cells)
for i, cell in enumerate(cells):
_, side, col, row = cell
for face_id in range(9):
n_addr = get_neighbor_for_face(side, col, row, face_id)
if n_addr in all_cells_set:
Topology[i, face_id] = cell_to_index[n_addr]
_, n_side, n_col, n_row = n_addr
for rev_f in range(9):
if get_neighbor_for_face(n_side, n_col, n_row, rev_f) == (0, side, col, row):
Neighbor_Faces[i, face_id] = rev_f
break
# Инициализация состояний
States = np.full(N_CELLS, CELL_EMPTY, dtype=np.int32)
# В оригинале было 4 пчелы, но для теста оставлена 1 центральная.
# Если нужно 4, раскомментируйте строку ниже.
# bee_cells = [(0, 'Y', 0, 0)]
bee_cells = [(0, 'Y', 0, 0), (0, 'Y', 0, 1), (0, 'Y', 0, 2), (0, 'Y', 2, -1)]
bee_idx = [cell_to_index[cell] for cell in bee_cells]
for idx in bee_idx:
States[idx] = CELL_BEE
# ==========================================================
# 3. РАСЧЕТ НАЧАЛЬНЫХ ТЕМПЕРАТУР (Тепловой баланс отдельной пчелы)
# ==========================================================
def metabolic_power(T):
return max((29.133 - 0.739 * T) * 1e-3, Pmin)
def calculate_initial_temperatures(T_air_val):
l_val = (1.3493 * np.exp(-0.058 * T_air_val)) * 1e-3
R_th_hd_val = l_val / (k_hd * s_r)
R_th_ab_val = l_val / (k_ab * s_r)
def equations(vars_val):
Thd_v, Tth_v, Tab_v = vars_val
q_evap = ke0 * Tth_v
P_rest = metabolic_power(T_air_val)
p_s_hd = -eps * sigma * ((Thd_v - t0)**4 - (T_air_val - t0)**4) * A_hd
p_s_th = -eps * sigma * ((Tth_v - t0)**4 - (T_air_val - t0)**4) * A_th
p_s_ab = -eps * sigma * ((Tab_v - t0)**4 - (T_air_val - t0)**4) * A_ab
eq1 = p_s_hd + (Tth_v - Thd_v)/R_th_hd_val - (Thd_v - T_air_val)/R_head_air_open - q_evap
eq2 = P_rest + p_s_th - (Tth_v - Thd_v)/R_th_hd_val - (Tth_v - Tab_v)/R_th_ab_val - (Tth_v - T_air_val)/R_thorax_air_open
eq3 = p_s_ab + (Tth_v - Tab_v)/R_th_ab_val - (Tab_v - T_air_val)/R_ab_air_open - q_evap
return [eq1, eq2, eq3]
return fsolve(equations, [T_air_val + 10, T_air_val + 15, T_air_val + 10], xtol=1e-8)
Thd_init, Tth_init, Tab_init = calculate_initial_temperatures(T_ambient)
T_bee_avg_init = (m_head * Thd_init + m_thorax * Tth_init + m_ab * Tab_init) / (m_head + m_thorax + m_ab)
# ==========================================================
# 4. ПОДГОТОВКА МАТРИЦ ДЛЯ NUMBA
# ==========================================================
C_grid = np.zeros((N_CELLS, 13), dtype=np.float64)
for i in range(N_CELLS):
C_grid[i, 4:13] = C_wall_base
if States[i] == CELL_BEE:
C_grid[i, I_HD] = C_hd
C_grid[i, I_TH] = C_th
C_grid[i, I_AB] = C_ab
C_grid[i, I_AI] = C_air_bee
else:
C_grid[i, I_AI] = C_air_cell
# Двойной буфер для температур (избегаем аллокаций внутри цикла!)
T_grid_A = np.full((N_CELLS, 13), T_ambient, dtype=np.float64)
T_grid_B = np.full((N_CELLS, 13), T_ambient, dtype=np.float64)
for t_idx in bee_idx:
T_grid_A[t_idx, I_HD] = Thd_init
T_grid_A[t_idx, I_TH] = Tth_init
T_grid_A[t_idx, I_AB] = Tab_init
T_grid_B[t_idx, I_HD] = Thd_init
T_grid_B[t_idx, I_TH] = Tth_init
T_grid_B[t_idx, I_AB] = Tab_init
# ==========================================================
# 5. ЯДРО СИМУЛЯЦИИ (NUMBA JIT)
# ==========================================================
# БЫЛО:
# @njit(fastmath=True, cache=True)
# def do_cluster_step(T_cur, T_next, Topology, Neighbor_Faces, States, C_grid, dt, T_ambient):
# СТАЛО:
@njit(fastmath=True, cache=True)
def do_cluster_step(T_cur, T_next, Topology, Neighbor_Faces, States, C_grid, dt, T_ambient, heater_idx):
Q_entrance_Y = 0.0
Q_entrance_X = 0.0
half_cells = len(T_cur) // 2
for i in range(len(T_cur)):
state = States[i]
Thd = T_cur[i, I_HD]
Tth = T_cur[i, I_TH]
Tab = T_cur[i, I_AB]
# Расчет температуры воздуха в ячейке (Tai)
if state == CELL_BEE:
num = Tth/R_th_ai + Thd/R_hd_ai + Tab/R_ab_ai + T_ambient/R_ab_
den = 1.0/R_th_ai + 1.0/R_hd_ai + 1.0/R_ab_ai + 1.0/R_ab_
else:
R_ent = 1.0 / (h_entrance * A_entrance)
num = T_ambient / R_ent
den = 1.0 / R_ent
for f in range(9):
num += T_cur[i, 4 + f] / R_ai_wall_arr[f]
den += 1.0 / R_ai_wall_arr[f]
Tai = num / den
T_next[i, I_AI] = Tai
is_side_Y = (i < half_cells)
# Расчет теплопотерь
if state == CELL_BEE:
Q_cell = hA_ * (Tai - T_ambient)
Q_ab = (Tab - T_ambient) / R_ab_
if is_side_Y: Q_entrance_Y += (Q_cell + Q_ab)
else: Q_entrance_X += (Q_cell + Q_ab)
else:
Q_cell = hA * (Tai - T_ambient)
if is_side_Y: Q_entrance_Y += Q_cell
else: Q_entrance_X += Q_cell
# Расчет температур стенок
T_wall_avg = 0.0
for f in range(9):
T_wall_avg += T_cur[i, 4 + f]
T_wall_avg /= 9.0
p_s_th = -se_th * ((Tth - t0)**4 - (T_wall_avg - t0)**4)
p_s_hd = -se_hd * ((Thd - t0)**4 - (T_wall_avg - t0)**4)
p_s_ab = -se_ab_ * ((Tab - t0)**4 - (T_wall_avg - t0)**4)
p_s_total = p_s_th + p_s_hd + p_s_ab
for f in range(9):
Q_from_air = (Tai - T_cur[i, 4 + f]) / R_ai_wall_arr[f]
neighbor_idx = Topology[i, f]
if neighbor_idx != -1:
n_face = Neighbor_Faces[i, f]
Q_from_neighbor = (T_cur[neighbor_idx, 4 + n_face] - T_cur[i, 4 + f]) / R_wall_self[f]
else:
Q_from_neighbor = (T_ambient - T_cur[i, 4 + f]) / R_wall_self[f]
p_s_to_wall_k = ((A_lateral if f < 6 else A_bottom) / A_wall_total) * p_s_total if state == CELL_BEE else 0.0
T_next[i, 4 + f] = T_cur[i, 4 + f] + (dt / C_grid[i, 4 + f]) * (Q_from_air + Q_from_neighbor + p_s_to_wall_k)
# Расчет температур органов пчелы
if state == CELL_BEE:
P_rest = (29.133 - 0.739 * Tai) * 1e-3
if P_rest < Pmin: P_rest = Pmin
# Добавляем мощность "печки" только для центральной пчелы (индекс 0 в bee_idx)
# БЫЛО: # P_local = P_warm_up if i == bee_idx[0] else 0.0
# СТАЛО: P_local = P_warm_up if i == heater_idx else 0.0
l_dyn = (1.3493 * np.exp(-0.058 * Tai)) * 1e-3
R_th_hd_dyn = l_dyn / (k_hd * s_r)
R_th_ab_dyn = l_dyn / (k_ab * s_r)
T_next[i, I_TH] = Tth + (dt / C_grid[i, I_TH]) * (
P_rest + P_local + p_s_th - (Tth - Thd)/R_th_hd_dyn - (Tth - Tab)/R_th_ab_dyn - (Tth - Tai)/R_th_ai)
T_next[i, I_HD] = Thd + (dt / C_grid[i, I_HD]) * (
(Tth - Thd)/R_th_hd_dyn - (Thd - Tai)/R_hd_ai - ke0 * Tth)
T_next[i, I_AB] = Tab + (dt / C_grid[i, I_AB]) * (
p_s_ab + (Tth - Tab)/R_th_ab_dyn - (Tab - Tai)/R_ab_ai - (Tab - T_ambient)/R_ab_ - ke0 * Tth)
return Q_entrance_Y, Q_entrance_X
# ==========================================================
# 6. ЦИКЛ СИМУЛЯЦИИ И СБОР ДАННЫХ
# ==========================================================
dt_sim = 0.1
time_total = 600.0
steps = int(time_total / dt_sim)
# Предаллокация списков для логов
log_time = []
log_T_hd, log_T_center, log_T_ab, log_T_ai_bee, log_T_wall_bee, log_T_bee_avg = [], [], [], [], [], []
log_Power, log_Q_Y, log_Q_X = [], [], []
print(f"Задана сетка {cols}x{rows}. Запуск скоростной симуляции...")
t_idx = bee_idx[0]
# Прогрев Numba (первый вызов всегда медленный из-за компиляции)
do_cluster_step(T_grid_A, T_grid_B, Topology, Neighbor_Faces, States, C_grid, dt_sim, T_ambient, bee_idx[0])
start_time = time.time()
T_cur, T_next = T_grid_A, T_grid_B
for step in range(steps):
# Считаем шаг, записывая результаты из T_cur в T_next
Q_Y, Q_X = do_cluster_step(T_cur, T_next, Topology, Neighbor_Faces, States, C_grid, dt_sim, T_ambient, bee_idx[0])
# Меняем буферы местами (ссылки, а не копирование данных!)
T_cur, T_next = T_next, T_cur
if step % 10 == 0:
t = step * dt_sim
log_time.append(t)
log_T_hd.append(T_cur[t_idx, I_HD])
log_T_center.append(T_cur[t_idx, I_TH])
log_T_ab.append(T_cur[t_idx, I_AB])
log_T_ai_bee.append(T_cur[t_idx, I_AI])
log_T_wall_bee.append(T_cur[t_idx, 4])
T_bee_avg = (m_head * T_cur[t_idx, I_HD] + m_thorax * T_cur[t_idx, I_TH] + m_ab * T_cur[t_idx, I_AB]) / (m_head + m_thorax + m_ab)
log_T_bee_avg.append(T_bee_avg)
P_rest = max((29.133 - 0.739 * T_cur[t_idx, I_AI]) * 1e-3, Pmin)
log_Power.append((P_rest + P_warm_up) * 1000)
log_Q_Y.append(Q_Y * 1000)
log_Q_X.append(Q_X * 1000)
time_calc = time.time() - start_time
print(f"Общее время: {time_calc:.3f} c ({steps/time_calc:.0f} шагов/сек)")
# ==========================================================
# 7. ВИЗУАЛИЗАЦИЯ
# ==========================================================
# График 1: Динамика температур
plt.figure(figsize=(8, 5))
plt.scatter(log_time, log_T_center, s=2, c='r', label='Торакс', zorder=5)
plt.plot(log_time, log_T_hd, 'b-', label='Голова', linewidth=1.5)
plt.plot(log_time, log_T_ab, 'g-', label='Брюшко', linewidth=1.5)
plt.plot(log_time, log_T_ai_bee, 'c--', label='Воздух в ячейке', linewidth=1.5)
plt.plot(log_time, log_T_wall_bee, 'm--', label='Стенка ячейки', linewidth=1.5)
plt.plot(log_time, log_T_bee_avg, 'k:', label='Средняя T пчелы', linewidth=2)
plt.axhline(T_ambient, color='gray', linestyle=':', label=f'Наружный воздух ({T_ambient}°C)')
plt.xlabel('Время (с)'); plt.ylabel('Температура (°C)')
plt.title(f'Тепловой режим пчелы в кластере. Вариант {variant}')
plt.legend(); plt.grid(True, linestyle=':'); plt.tight_layout(); plt.show()
# График 2: Баланс мощности
fig, ax1 = plt.subplots(figsize=(8, 5))
ax1.plot(log_time, log_Power, 'orange', linewidth=2, label='Мощность термогенеза')
ax1.set_xlabel('Время (с)'); ax1.set_ylabel('Мощность (мВт)', color='orange')
ax2 = ax1.twinx()
ax2.plot(log_time, log_Q_Y, 'steelblue', linewidth=2, label='Потери в улочку Y')
ax2.plot(log_time, log_Q_X, 'crimson', linewidth=2, label='Потери в улочку ^')
ax2.set_ylabel('Теплопотери (мВт)', color='k')
lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper right')
ax1.grid(True, linestyle=':'); ax1.set_title('Баланс мощности: генерация vs потери'); plt.show()
# График 3: 2D-карта температур сот
temp_Y_clean, temp_X_clean = {}, {}
for i, cell in enumerate(cells):
_, side, col, row = cell
T_air_final = T_cur[i, I_AI]
if side == 'Y': temp_Y_clean[(col, row)] = T_air_final
else: temp_X_clean[(col, row)] = T_air_final
size = 1.0
width = np.sqrt(3) * size
dy = 1.5 * size
def hex_corner(center, i):
angle_rad = np.pi / 180 * (60 * i - 30)
return (center[0] + size * np.cos(angle_rad), center[1] + size * np.sin(angle_rad))
def primal_to_pixel(col, row): return (width * (col + 0.5 * (row & 1)), -dy * row)
def dual_to_pixel(col, row): return (width * (col + 0.5 * (1 - (row & 1))), -dy * row - size / 2)
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_aspect('equal')
all_temps = list(temp_Y_clean.values()) + list(temp_X_clean.values())
vmin, vmax = min(all_temps), max(all_temps)
for col in range(-int(cols/2), int(cols/2)+1):
for row in range(-int(rows/2), int(rows/2)+1):
if (col, row) in temp_Y_clean:
x, y = primal_to_pixel(col, row)
corners = [hex_corner((x, y), i) for i in range(6)]
hex_x, hex_y = zip(*corners)
color = plt.cm.viridis((temp_Y_clean[(col, row)] - vmin) / max(vmax - vmin, 1e-5))
ax.fill(hex_x, hex_y, color=color, edgecolor='steelblue', linewidth=1.2, alpha=0.8)
if (col, row) == (0, 0): ax.fill(hex_x, hex_y, facecolor='none', edgecolor='black', linewidth=3)
ax.text(x, y, f"{temp_Y_clean[(col, row)]:.2f}", ha='center', va='center', fontsize=7, color='b')
for col in range(-int(cols/2), int(cols/2)+1):
for row in range(-int(rows/2), int(rows/2)+1):
if (col, row) in temp_X_clean:
x, y = dual_to_pixel(col, row)
corners = [hex_corner((x, y), i) for i in range(6)]
hex_x, hex_y = zip(*corners)
ax.plot(hex_x + (hex_x[0],), hex_y + (hex_y[0],), color='crimson', linewidth=1.2)
ax.text(x, y, f"{temp_X_clean[(col, row)]:.2f}", ha='center', va='center', fontsize=7, color='crimson')
ax.axis('off'); ax.autoscale_view()
sm = plt.cm.ScalarMappable(cmap=plt.cm.viridis, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm.set_array([])
plt.colorbar(sm, ax=ax, shrink=0.8).set_label('Температура воздуха (°C)', fontsize=10)
plt.title(f'Температурное поле сот. Вариант {variant}\n(Синяя — сторона "Y", Красная — сторона "^")')
plt.tight_layout(); plt.show()
# Финальный вывод в консоль
print(f"\nФинальные результаты для центральной пчелы:")
print(f"Торакс: \t\t{log_T_center[-1]:.2f}°C (было {Tth_init:.2f}°C)")
print(f"Средняя пчела: \t{log_T_bee_avg[-1]:.2f}°C (было {T_bee_avg_init:.2f}°C)")
print(f"Воздух в ячейке: \t{log_T_ai_bee[-1]:.2f}°C")
print(f"Мощность: \t\t{log_Power[-1]:.2f} мВт")
Комментарии
Отправить комментарий