Клеточный автомат тепла и массы — 2D-разрез хребта (прототип)

Проверка идеи пользователя (карточка docs/archive/plan/heat-ca-prototype.md): вместо поля ветра «сразу» — клетки, которые обмениваются теплом и массой с соседями, и смотрим, что куда течёт над прогретым склоном. Игру не трогает. Физика и планы — docs/research/slope_wind.md (§8–§11), docs/archive/plan/wind-field.md.

Запуск

cd tools/research/heat_ca
uv venv .venv && uv pip install --python .venv/bin/python -r requirements.txt   # один раз
                                                    # (python3 -m venv без ensurepip здесь не работает)
.venv/bin/python run.py s1_sun_one_slope            # один сценарий (клетка 50 м, GPU)
.venv/bin/python run.py all                         # все сценарии
.venv/bin/python run.py all --cell 100              # клетка — параметр: 200/100/50/25/12.5/6.25 м
.venv/bin/python run.py s1_sun_one_slope --pressure acoustic   # чисто локальный режим (см. ниже)
.venv/bin/python run.py s1_sun_one_slope --device cpu          # без GPU (numpy, для отладки)
.venv/bin/python study.py cells [--pressure acoustic]  # сходимость по клетке + время на GPU + оценка 3D
.venv/bin/python study.py solvers                      # сколько итераций нужно выравниванию массы

Сценарии: s1_sun_one_slope (утреннее солнце 25° слева на крутой склон, пологий — в косом свете), s2_both_slopes (оба склона прогреты одинаково), s2_sym (то же на симметричном хребте — проверка симметрии схемы), s3_wind (как s2 + фоновый ветер 3 м/с слева направо), s4_inversion (как s2 + инверсия 1500–1800 м над подножием). Хребет 500 м, крутой склон слева (до ~33°), пологий справа (до ~16°). Разрез 6,4 × 3,2 км (с ветром — 9,6 км: запас по ветру).

GPU — CuPy (CUDA) только для прототипа: в игре такой счёт пойдёт вычислительными шейдерами Vulkan (нужна работа на AMD), см. docs/archive/plan/wind-field.md. --device cpu оставлен для отладки, время на CPU не сравниваем.

Что в клетке и какие правила (простыми словами)

В клетке — масса m (норма 1) и тепло H = m·θ (θ — потенциальная температура: температура, приведённая к одной высоте, чтобы «теплее соседа сверху» сравнивалось честно). Между каждой парой соседних клеток — поток через их общую грань: сколько ушло из одной, столько пришло в другую. Поэтому масса и тепло сохраняются точно (до округления), а скорость воздуха — это и есть поток через грани (для стрелок — среднее по граням клетки).

За шаг времени:

1. Нагрев. Земля отдаёт тепло в нижнюю воздушную клетку столбца (поток ∝ солнцу на склон).
2. Тёплое поднимается. Поток вверх через грань ускоряется, если клетки у грани теплее фона на своей высоте (фон — температура окружающего воздуха, падающая с высотой). Никакого фиксированного «вверх больше»: поднятый воздух охлаждается (несёт свою θ), и где фон теплеет с высотой (инверсия), он быстро оказывается холоднее соседей — подъём гаснет сам.
3. Инерция и трение. Потоки переносят сами себя, вязкость сглаживает их между соседями, у земли — трение; у потолка — «губка», гасящая волны.
4. Недостача массы подсасывает соседей. Если потоки унесли из клетки больше, чем принесли, в ней «давление» падает, и потоки с соседями поправляются в её сторону. Два режима:
   • mg (по умолчанию): поправка сразу такая, чтобы недостачи не осталось, — это уравнение Пуассона для давления; итерации «каждая клетка смотрит на соседей» (Якоби) ускорены многосеточным методом (1–2 V-цикла на шаг, с тёплым стартом);
   • acoustic: буквально идея пользователя, без глобального решателя — избыток массы клетки сам даёт «давление» c²·(m−1), потоки разгоняются перепадом между соседями, масса переносится этими потоками; несколько подшагов на шаг («медленный звук» c = 60 м/с ≫ ветра) и гашение дивергенции. Недостача массы в клетке остаётся ~0,1–2 % (∝ 1/c²), картина та же.
5. Перенос. Масса и тепло переходят через грани по потокам (тепло — с температурой клетки, откуда дует), плюс теплопроводность между соседями по отклонению от фона, плюс слабое выхолаживание к фону (τ = 2 ч — явный сток, чтобы было установившееся состояние; в балансе).

Это и есть конвекция Буссинеска на сетке «клетки + грани» (MAC), записанная локально: п. 2 — плавучесть относительно фона, п. 4 — давление, итерации п. 4 — Якоби для уравнения Пуассона.

Шаг по времени — устойчивый: dt ≤ 0,35·клетка/8 м/с (перенос) и ≤ 0,2·клетка²/K (теплопроводность); «пилы» нет. Установление — когда за 10 мин модельного времени скорость меняется < 0,02 м/с и температура < 0,03 К.

Границы: в штиль бока — стенки (замкнутая долина между хребтами: что поднялось, опускается внутри разреза); с ветром — слева заданный приток, справа выход «со сносом», по бокам губки. У краёв разреза нагрев сходит на нет (край — не склон). С ветром трение действует на отклонение от фонового профиля (фон считается уравновешенным крупным перепадом давления, иначе трение тормозит набегающий поток ещё в долине и поднимает его — ложный подъём).

Картинки (out/<сценарий>/)

• temp.png — θ′: насколько воздух теплее (+) / холоднее (−) фона на той же высоте; рельеф — коричневый, контур — «лесенка» клеток модели;
• flow.png — стрелки ветра через 200 м (длина — скорость, ключ внизу), цвет стрелки — вертикальная скорость (красное — подъём, синее — опускание), серое под стрелками — прогретый воздух;
• w_profile.png — w по X на 50/150/300/600 м над рельефом + рельеф + поток нагрева;
• mass_balance.png — масса и тепло: изменение = притоки через границы + нагрев + стоки; невязки (лог-шкала) против допуска 0,1 %;
• evolution.mp4 — как складывается циркуляция от старта до установления;
• s3: flow_no_heating.png (тот же ветер без нагрева — чистое обтекание) и flow_thermal.png (разность: что добавил нагрев);
• summary.md — числа и условия приёмки; result.json — то же для машин.
• out/convergence*/ — сходимость по клетке (conv_w_profile.png, conv_flow.png); out/solvers/solvers_w300.png — способы выравнивания массы; out/summary.md — общий итог.

Итоги

Ниже — копия out/summary.md (общий вывод, таблицы времени и сходимости).

Прототип по карточке docs/archive/plan/heat-ca-prototype.md. Код, запуск и смысл картинок — README.md. Хребет 500 м, крутой склон слева (до ~33°), пологий справа (до ~16°); разрез 6,4 × 3,2 км; фон почти безразличный (θ растёт на 1 К/км, т. е. T падает на ~8,8 К/км); нагрев склонов до 250 Вт/м². Счёт — GPU (RTX 4070 SUPER, CuPy + CUDA Graph, float32); в игре это будут вычислительные шейдеры Vulkan (должно работать на AMD), CuPy — только для прототипа.

Правила автомата (простыми словами)

1. Земля греет нижнюю клетку столбца (по солнцу на склон).
2. Поток вверх через грань ускоряется, если клетки теплее фона на своей высоте; поднятый воздух несёт свою температуру — на инверсии быстро становится холоднее соседей, подъём гаснет сам.
3. Потоки переносят сами себя, вязкость и трение сглаживают; губка у потолка.
4. Недостача массы в клетке «подсасывает» соседей (давление) — многосеточно (mg) или чисто локально (acoustic, «медленный звук», без глобального решателя).
5. Масса и тепло переходят через грани парными потоками (что ушло из одной — пришло в другую) → сохранение точное; плюс теплопроводность, плюс слабое выхолаживание к фону (τ = 2 ч, явный сток).

Это конвекция Буссинеска, записанная локально: п. 2 — плавучесть относительно фона, п. 4 — уравнение Пуассона для давления (итерации «клетка смотрит на соседей» = Якоби).

Сценарии: результат против условий приёмки (клетка 50 м)

• Сохранение (п. 1): невязка массы ~10⁻¹¹ от массы области, тепла ~10⁻⁵ от суммарного нагрева (float32); допуск 0,1 % выполнен с запасом на порядки. Недовыравненная масса в клетке ~2·10⁻⁸ (многосеточно) или 0,1–2 % (локальный режим, ∝ 1/c²). Графики — mass_balance.png.
• П. 2 (сценарий 1): подъём над прогретым склоном и гребнем, приток у подножия к склону, возвратный поток на 0,8–1,3 км и опускание над долиной — s1_sun_one_slope/flow.png, evolution.mp4.
• П. 3: схождение и столб подъёма над гребнем; на несимметричном хребте столб сдвинут на ~0,5 км к пологому склону (с него приходит больше тёплого воздуха); на симметричном (s2_sym) — строго над гребнем, картина зеркальная до 10⁻⁵.
• П. 4 (инверсия): подъём гаснет у нижней кромки инверсии (1550 м против 1800 без неё), под инверсией — растекание в стороны до 1,6 м/с и холодная «шапка» перебега; выше инверсии |w| ≤ 0,03 м/с, |θ′| ≤ 0,06 К — s4_inversion/flow.png, temp.png.
• П. 5 (ветер): картина не разваливается; полное поле — это прежде всего обтекание хребта устойчивым воздухом (горная волна: подъём на наветренном склоне и над долиной перед ним, сильное опускание за гребнем). Тепловая часть (с нагревом минус без) — наклонённый по ветру подъём от наветренного склона через гребень (dx/dz ≈ +0,5), тёплый слой сдувается за хребет — s3_wind/flow_thermal.png, flow_no_heating.png. При 4 м/с тепловую часть почти не видно — её забивает волна (честно: при ветре 4+ м/с и такой устойчивости тепло склона — поправка ~0,3 м/с).
• П. 6: «пилы» нет, шаг устойчивый (перенос: dt ≤ 0,35·клетка/8 м/с; теплопроводность: ≤ 0,2·клетка²/K), установление за ~200 мин модельного времени — 1,2–15 тыс. шагов по клетке.

Похоже ли на слова пилота. Да по всем трём пунктам карточки: подъём у прогретого склона и над гребнем (максимум — над верхней частью склона, 0,4–0,7 км над землёй), приток снизу к подножию (~0,5 м/с, вверх по склону ~0,5–0,8 м/с — склоновый ветер), потолок на инверсии (столб останавливается у её нижней кромки и растекается). Опускание над долиной слабое и широкое (−0,1…−0,2 м/с), рядом со столбом — сильнее (до −0,75 м/с). Числа порядка настоящих склоновых циркуляций; сила зависит от вязкости/теплопроводности (30 м²/с — грубое «перемешивание»).

Сходимость по клетке и время (сценарий 1)

Сходимость по клетке и время на GPU (NVIDIA GeForce RTX 4070 SUPER, CuPy + CUDA Graph, float32; сценарий 1; выравнивание массы — многосеточный Пуассон, 2 V-цикла на шаг)

Оценка для 3D (квадрат 40×40 км × 32 уровня, та же физика, тот же GPU): мс/шаг = max(0.49 мс — пол накладных расходов, 2.88 нс на клетку·шаг по самой большой 2D-сетке × 1.5 за 3D); шагов — как в 2D при той же клетке; память — 432 байт на клетку × 1,3.

Прототип держит в пуле CuPy все временные массивы шага (CUDA Graph) — отсюда ~430 байт на клетку; в вычислительном шейдере нужно ~12–16 чисел float32 на клетку (u, v, w, m, H, p, правая часть, невязка, уровни многосеточного) ≈ 64 байт.

Сходимость по клетке и время на GPU (NVIDIA GeForce RTX 4070 SUPER, CuPy + CUDA Graph, float32; сценарий 1; выравнивание массы — чисто локально («медленный звук» 60 м/с))

Оценка для 3D (квадрат 40×40 км × 32 уровня, та же физика, тот же GPU): мс/шаг = max(0.35 мс — пол накладных расходов, 2.73 нс на клетку·шаг по самой большой 2D-сетке × 1.5 за 3D); шагов — как в 2D при той же клетке; память — 432 байт на клетку × 1,3.

Прототип держит в пуле CuPy все временные массивы шага (CUDA Graph) — отсюда ~430 байт на клетку; в вычислительном шейдере нужно ~12–16 чисел float32 на клетку (u, v, w, m, H, p, правая часть, невязка, уровни многосеточного) ≈ 64 байт.

Сходимость по клетке (out/convergence/conv_w_profile.png, conv_flow.png): картина (где подъём, приток, возвратный поток, опускание, высота подъёма 1400 м) одинакова от 200 до 6,25 м. Меняются:

• макс. подъём растёт с измельчением (1,21 → 1,42 → 1,65 → 1,81 → 1,90 → 1,95 м/с) — сходится как первый порядок, 50 м даёт ~85 % от предела;
• склоновый ветер вверх по склону: на 200 м его нет (+0,02 м/с — склон 2 клетками), на 100 м уже +0,47 (предел ~0,67);
• приток у подножия на грубых сетках завышен (0,71 на 200 м против ~0,40). Рекомендация: самая грубая клетка с правильной картиной — 100 м (есть склоновый ветер, подъём в нужном месте, ~75 % силы); 200 м годится только как фон (циркуляция долины без склонового слоя). Для 3D при загрузке — 200 м по всей области и 50–100 м в окне у пилота (как WF-14).

Время. Число шагов до установления ∝ 1/клетка (шаг по времени ∝ клетке, а устанавливается ~3 часа модельного времени — медленный прогрев долины). На малых сетках GPU ждёт запуски (пол ~0,4–0,5 мс/шаг даже с CUDA Graph), на 3,125 м (2 млн клеток) — 2,9 нс на клетку·шаг. Оценка 3D 40×40 км × 32 уровня: 200 м — ~7 с, 100 м — ~1 мин, 50 м — ~8 мин, 25 м — ~1 ч (тот же GPU; на слабой AMD — в 3–10 раз дольше). Память в компактной реализации ~64 байт/клетку: 78 МБ (200 м) … 1,25 ГБ (50 м).

Выравнивание массы: сколько итераций нужно

Выравнивание массы: сколько итераций нужно (сценарий 1, клетка 50 м, GPU)

| Способ | Итераций (подшагов) на шаг | Макс. |m−1| в клетке (недовыравнено) | Шагов до установления | мс/шаг | Макс. подъём, м/с | Приток у подножия, м/с | Высота подъёма, м | |—|—|—|—|—|—|—|—| | многосеточный V-цикл (Пуассон) | 1 | 5.9e-08 | 5754 | 0.392 | 1.65 | +0.49 | 1400 | | многосеточный V-цикл (Пуассон) | 2 | 2.2e-08 | 5754 | 0.530 | 1.65 | +0.49 | 1400 | | многосеточный V-цикл (Пуассон) | 4 | 2.0e-08 | 5754 | 0.812 | 1.65 | +0.49 | 1400 | | Якоби, локальные итерации давления, снимать накопленный избыток массы | 50 | расходимость: |v|max=nan на шаге 274 | | | | | | | Якоби, локальные итерации давления, без возврата накопленной массы | 50 | 3.0e-03 | 5754 | 0.297 | 1.65 | +0.49 | 1400 | | Якоби, локальные итерации давления, без возврата накопленной массы | 500 | 8.0e-04 | 5754 | 0.715 | 1.65 | +0.49 | 1400 | | локально: «медленный звук» c = 30 м/с | 3 | 2.3e-02 | 5754 | 0.288 | 1.65 | +0.50 | 1400 | | локально: «медленный звук» c = 60 м/с | 6 | 5.6e-03 | 5754 | 0.355 | 1.64 | +0.51 | 1400 | | локально: «медленный звук» c = 120 м/с | 11 | 1.3e-03 | 5754 | 0.468 | 1.63 | +0.50 | 1400 |

• Многосеточному хватает 1 V-цикла на шаг с тёплым стартом (недостача массы 6·10⁻⁸).
• Буквальная «клетка смотрит на соседей» (Якоби), если каждый шаг пытаться вернуть накопленную недостачу массы, разваливается (двойной счёт недостачи → раскачка); без возврата работает, но масса в клетках «плывёт» (0,1–0,3 %), и на больших сетках Якоби сходится всё медленнее (∝ N²).
• Лучший локальный вариант — «медленный звук» (--pressure acoustic): недостача массы сама даёт давление, 3–11 подшагов по соседям на шаг, без глобального решателя; картина та же (подъём 1,63–1,65 м/с), быстрее многосеточного (0,29–0,47 против 0,39–0,53 мс/шаг), масса в клетке 0,1–2 %. Для вычислительного шейдера это самое удобное: только соседи, никаких уровней и редукций; на AMD переносимо.

Сравнение с решателем Пуассона из плана поля ветра

Выравнивание массы здесь — это то же уравнение Пуассона, что в docs/archive/plan/wind-field.md (∇·(∇λ) = −∇·u₀, те же проводимости граней, многосеточный метод). Разница в том, что поле ветра решает его один раз (диагностическое поле без инерции и плавучести: 15–30 итераций PCG), а автомат — на каждом шаге (1 V-цикл × тысячи шагов) плюс перенос тепла и импульса. То есть автомат = решатель поля ветра + плавучесть + время. Отсюда цена: при 200 м — секунды против долей секунды; при 50 м — минуты против секунды. Зато автомат даёт то, чего поле по массе не даёт в принципе (§11.2 исследования): приток к прогретому склону, столб над гребнем, возвратный поток и опускание, потолок на инверсии, наклон по ветру.

Вывод для игры

• Идея рабочая и физически честная: локальные правила дают правильную качественную картину и правильные порядки величин; сохранение точное; счёт устойчивый.
• Статичный прогон до установления при загрузке/смене часа реален только на грубой сетке: 3D 200 м — ~5–10 с на 4070 SUPER (на слабой AMD — до минуты), 100 м — ~1 мин (слишком долго для загрузки, можно фоном). 50 м по всей области — нет; 50 м в окне 8×8 км у пилота — ~30 с фоном.
• Ускорители: локальный «медленный звук» вместо многосеточного (проще шейдер, быстрее); старт не с покоя, а с прошлого часа (установление не с нуля, а на ~30–60 мин модельного времени — в 3–5 раз меньше шагов); грубая сетка для фона + мелкая в окне (WF-14).
• Что это даст пилоту поверх поля ветра по массе: склоновый ветер и приток снизу (трава/пыль тянутся к склону), широкое опускание над долиной, высоту «потолка» по инверсии, снос тепла ветром. Термики как пузыри и «провалы» — по-прежнему отдельная модель (автомат на 50–200 м их не разрешает).
• Рекомендация: не заменять план поля ветра, а рассматривать автомат как этап 2–3 связанной системы (§11.3): тот же решатель + плавучесть от солнца на склонах, 200 м по области, раз в игровой час, фоном; начать с 2D→3D прототипа на Vulkan с локальным режимом давления.