# Этап 2: модуль Власова (частицы) для 3D-расчётов — проект

Дата: 30.09.2026. Контекст: этап 1 (план §6.12) показал, что скалярное поле как прокси пыли непригодно для
режима малых μ: предел Каупа M_max = 0.633/m превышает массу первого горизонта 1D-карты при любой доступной m,
а стоимость ∝ m (осцилляции поля задают шаг). Физически правильный источник — бесстолкновительные частицы
(уравнение Власова), как в 1D-коде V1. Этот документ фиксирует выбор платформы, уравнения, объём работ и тесты.

## 1. Платформа

| Вариант | Частицы | AH finder | AMR | Вывод |
|---|---|---|---|---|
| GRChombo + Chombo (форк GRChombo) | нет: ParticleTools отсутствует в форке и в публичном Chombo 3.2 (ChomboRelease3.2.10); пришлось бы писать контейнеры, перераспределение по боксам/уровням, депозит | есть (PETSc), проверен | Chombo, субциклирование | много инфраструктурной работы |
| GRTeclyn (AMReX) | AMReX ParticleContainer, `ParticleToMesh`/`MeshToParticle`, перераспределение и уровни — готово; в GRTeclyn уже есть `ParticleInterpolator` (Lagrange 4-го порядка + производные полей в точках частиц) | ветка `AHFinder` активна (03.09.2026), `Source/AHFinder` | AMReX, субциклирование | **выбор**: минимум своей инфраструктуры |

Выбор: **GRTeclyn, ветка AHFinder** (или main после слияния). Риск: ветка не доведена — тогда AH для e = 0
считать по сферическим срезам самостоятельно (поиск R = 2M по лучам), для e ≠ 0 — ждать слияния.

## 2. Уравнения

Частица p с массой m_p (комовингный вес), положение x^i_p, ковариантный импульс u_i (пространственная
часть 4-скорости, нижние индексы), Γ = √(1 + γ^{ij}u_i u_j) (γ^{ij} = χ h̃^{ij}). Уравнения геодезических в 3+1
(Vincent+2012, как в 1D-коде V1, обобщение на 3D):

- dx^i/dt = α γ^{ij} u_j / Γ − β^i
- du_i/dt = −Γ ∂_i α + u_j ∂_i β^j − (α/2Γ) u_j u_k ∂_i γ^{jk}

Источники в правой части CCZ4 (нормальный наблюдатель), депозит на сетку с весом 1/(√γ dV) (√γ = χ^{−3/2}):

- ρ = Σ m_p Γ_p W(x − x_p) / (√γ dV)
- S_i = Σ m_p u_{i,p} W / (√γ dV)
- S_{ij} = Σ m_p u_{i,p} u_{j,p} W / (Γ_p √γ dV),  S = γ^{ij} S_ij

Ядро W: CIC/TSC (AMReX `ParticleToMesh` с пользовательской функцией депозита); на границах уровней —
депозит на своём уровне + усреднение вниз (AMReX `average_down`) и заполнение призрачных ячеек, как для полей.

Класс материи `ParticleMatter` для `CCZ4RHSWithMatter`: e-m тензор читается из сеточных переменных
(ρ, S_i, S_ij — 10 компонент, «эволюционные» с нулевой правой частью, чтобы интерполяция по уровням и
рестарт работали как для полей); частицы — отдельный `ParticleContainer` с толчком после каждого шага RK4
(или на каждой стадии — как в V1: KDK внутри стадий; для начала — толчок на шаге с полями,
взятыми в середине шага, что даёт 2-й порядок по времени для частиц; поля — 4-й).

## 3. Объём работ (оценка)

1. Сборка GRTeclyn (AHFinder) + AMReX на сервере, прогон примеров ScalarField/BinaryBH — 1 день.
2. `ParticleMatter`: депозит ρ, S_i, S_ij; e-m тензор; связи — 2 дня.
3. Толчок частиц: интерполяция α, β, χ, h̃, их производных в точках частиц (расширение
   `ParticleInterpolator`), интегратор, перераспределение по боксам/уровням после толчка — 3 дня.
4. Начальные данные: FLRW-пыль (равномерные частицы, u_i = 0), LTB при t₀ (генератор
   `ltb_to_grchombo.py` → частицы с весами m'(r)dr и u_i = 0 в синхронном срезе: связи точны, как для скаляра) — 1 день.
5. Тесты (см. §4) — 3 дня. Итого ≈ 2 недели чистой работы.

## 4. Тесты

1. FLRW-пыль: ρ ∝ a⁻³, K = −3H, невязка связей от дискретности ∝ 1/√N_p.
2. LTB e = 0, μ = 0.3 и 0.1 (холодная пыль): ρ(0,t), R(t) оболочек, первый AH и кривая M_AH(t) против
   1D-кода (таблица в devlog от 30.09) — цель: совпадение времени достижения 0.1 и 0.3 M_H в пределах 5 %.
3. Тёплая пыль σ_ent: та же карта против 1D (σ = 0.016, 0.054).
4. Затем e ≠ 0 (§6.12, этап 1b): e = 0.13, 0.19 при μ = 0.04–0.06 из ICs Зельдовича/PM-3D при t₀
   (связи — GRTresna по депонированным ρ, S_i).

## 5. Что уже есть и переиспользуется

- Генератор LTB-среза в изотропных координатах (`ltb_to_grchombo.py`), калибровка K − f(x)K_far (§ devlog 29–30.09),
  вложенные сферы AMR с порогом плотности, анализ M_AH(t) (`pbh_analyse.py`).
- 1D-код V1 как эталон (карты холодная/тёплая, `runs/v1/ref3d`).

## 6. Реализация: раскладка кода (GRTeclyn, ветка AHFinder fd76177 + AMReX development)

Каталог `pbhgr/grteclyn/PBHVlasov/` (копируется в `GRTeclyn/Examples/PBHVlasov/`):

| Файл | Назначение |
|---|---|
| `StateVariables.hpp` | CCZ4 + 10 сеточных компонент материи `c_rho_p, c_S1..3, c_S11..33` (чётности как у ρ, S_i, S_ij) |
| `ParticleMatterVars.hpp` | `class ParticleMatterVars : CCZ4Vars` — доступ к ρ, S_i, S_ij ячейки |
| `ParticleMatter.hpp` | `matter_t` для `CCZ4RHSWithMatter`: `compute_emtensor` читает сеточные ρ, S_i, S_ij; `add_matter_rhs` = 0 |
| `VlasovParticles.hpp/.cpp` | `amrex::ParticleContainer<NStructReal=4 (m, u_1, u_2, u_3)>`; `deposit(lev, MultiFab&)` через `amrex::ParticleToMesh` (CIC, веса m Γ/(√γ dV), m u_i/(√γ dV), m u_i u_j/(Γ √γ dV)); `push(lev, dt, state, deriv_state)` через `amrex::MeshToParticle` (CIC-интерполяция α, β^i, χ, h̃_ij и их ∂_i, предварительно посчитанных на сетке в отдельный MultiFab), интегратор RK2/KDK; `Redistribute()` после толчка |
| `PBHVlasovLevel.hpp/.cpp` | `specific_post_init`: частицы из таблицы (FLRW/LTB), депозит; `specific_eval_rhs`: как ScalarField, но `CCZ4RHSWithMatter<ParticleMatter>`; `specific_advance`: толчок частиц + депозит для следующего шага; `specific_post_regrid`: `Redistribute` и депозит на новых уровнях; диагностика: ⟨K⟩, ρ_c, Ham/Mom, AH (ветка AHFinder) |
| `Main_PBHVlasov.cpp`, `GNUmakefile`, `params.txt` | `USE_PARTICLES = TRUE`, `Src/Particle/Make.package` |

Порядок: (1) скелет с нулём частиц компилируется и даёт вакуум/FLRW при ручном ρ; (2) депозит + FLRW-пыль
(u_i = 0, равномерные частицы) — тест a(t), K(t); (3) толчок — тест: свободная частица в FLRW (u_i = const,
dx/dt = α γ^{ij}u_j/Γ), затем LTB e = 0 при μ = 0.3 против 1D (таблица от 30.09).
Депозит только на уровне частицы; на грубые уровни — `average_down`; сглаживание ядра: CIC на всех уровнях,
при необходимости TSC. Число частиц: ≥ 8 на ячейку самого мелкого уровня в ядре (в 1D нужны были 32 оболочки
на ячейку для сходимости релаксации пересечений — здесь 3D-аналог проверяется сходимостью по N_p).