# 数値手法 ## 1. 数値単位 多くのシナリオは無次元のcode unitを用います。単位の対応はシナリオごとに異なり、現実の天体を定量予測することを目的としていません。 宇宙論モードはboxを`Mpc/h`、速度をpeculiar `km/s`、粒子質量を`Msun/h`として表示します。ただし初期条件と時間発展は教育用近似です。 ## 2. 重力多体系計算 ### 2.1 平面拘束または3次元Newton重力 粒子`i`の加速度を ```text a_i = G Σ_{j≠i} m_j (r_j-r_i) / (|r_j-r_i|²+ε²)^(3/2) ``` で計算します。`ε`はPlummer型softeningです。 2次元modeでも`newton-3d`を選んだ場合、座標を平面へ拘束するだけで力則は3次元Newton重力です。薄い円盤、惑星軌道、2D表示の星団に適します。 ### 2.2 真の2D Poisson点重力 2次元Poisson方程式に対応する対数potentialでは、点質量の力の大きさは概ね`1/r`です。 ```text Φ ∝ ln r |F| ∝ 1/r ``` 3次元Newton重力の軽量版ではありません。次元依存性を学ぶ教材用です。 ### 2.3 Direct summation 全pairを計算します。 ```text pair count = N(N-1)/2 complexity ≈ O(N²) ``` 小粒子数、collisional toy model、基準解に利用します。Worker版はpair-rangeを総pair数が均等になるよう分割し、部分加速度をcoordinatorで和します。 ### 2.4 Barnes–Hut 2Dではquadtree、3Dではoctreeを構築します。遠方nodeについて、node size `s`、距離`d`、opening angle `θ`を用いて一括近似します。 ```text s/d < θ ならnodeを一つの質点として近似 ``` - 小さい`θ`: 高精度・低速 - 大きい`θ`: 低精度・高速 公開実装はmonopole近似です。open boundaryのみを対象とし、periodic gravityではPMを使います。 ### 2.5 Particle–Mesh 周期boxでは次の順に計算します。 1. Cloud-in-Cellで粒子質量をmeshへdeposit 2. 平均密度を差し引いてdensity contrastを作る 3. radix-2 FFT 4. Fourier空間でPoisson方程式を解く 5. spectral derivativeでforce fieldを作る 6. inverse FFT 7. CICで粒子位置へ補間 3D cell数は`Nmesh³`です。mesh一辺を2倍にするとcell数は8倍になります。 PMは長距離周期重力と大粒子数に向きますが、cell幅より小さいforceは解像しません。本アプリはTreePMではありません。 ## 3. 時間積分 ### 3.1 Forward Euler ```text v_{n+1} = v_n + a_n dt x_{n+1} = x_n + v_n dt ``` 一次精度で、保存系の長時間積分ではenergy driftを生じやすい比較用手法です。 ### 3.2 Kick–Drift–Kick leapfrog ```text v_{n+1/2} = v_n + (dt/2) a(x_n) x_{n+1} = x_n + dt v_{n+1/2} v_{n+1} = v_{n+1/2} + (dt/2) a(x_{n+1}) ``` 二次精度、time reversible、symplecticな積分法です。エネルギー誤差が長時間で有界に振動しやすい特徴があります。 ### 3.3 Global adaptive timestep 最大加速度とsofteningからglobal dtを縮小します。全粒子が同じdtを使います。individual timestepではありません。 ## 4. 衝突・合体 `contactModel=merge`では、距離がcollision radius条件を満たした二粒子を完全非弾性合体します。 保存する量: - mass - linear momentum - center-of-mass position 一般に保存しない量: - kinetic energy 破砕、反発、材料強度、spinは含みません。 ## 5. 宇宙論初期条件 ### 5.1 粒子load 一辺`n`の3D Lagrangian latticeを置きます。 ```text Nparticle = n³ ``` `n`は8–96の整数です。FFT meshはradix-2に制限されますが、粒子格子は2のべき乗である必要はありません。 ### 5.2 Gaussian random field seeded real white noiseをFFTし、smooth matter transfer shape、primordial tilt `n_s`、WDM suppressionを適用します。 ### 5.3 Discrete sigma8 normalization box内の離散Fourier modeに対して、半径`8 Mpc/h`のtop-hat windowを使い、設定`σ8`へ規格化します。有限boxと離散modeのため、連続無限体積の厳密な`σ8`とは異なります。 ### 5.4 1LPT / Zel'dovich approximation 変位field`Ψ`を使い、Lagrangian coordinate `q`から ```text x(q,a) = q + D(a) Ψ(q) ``` として初期位置を作ります。速度は線形成長率の近似から与えます。 2LPTではないため、初期transientを精密に抑える用途には不十分です。 ### 5.5 Comoving evolution scale factorを時間変数として周期PM重力で進めます。背景膨張は選択した`Ωm, ΩΛ, Ωr, w0`から計算します。 ## 6. SPH ### 6.1 密度推定 ```text ρ_i = Σ_j m_j W(|r_i-r_j|, h) ``` 2Dと3Dではkernel normalizationが異なります。 - 2D: surface density - 3D: volume density ### 6.2 Cubic-spline kernel support radiusは`2h`です。自動試験では2D・3Dの数値積分でkernel normalizationを確認します。 ### 6.3 Pressure force 対称形のpressure gradientを用います。 ```text dv_i/dt = -Σ_j m_j (P_i/ρ_i² + P_j/ρ_j² + Π_ij) ∇W_ij ``` `Π_ij`は人工粘性です。 ### 6.4 Monaghan型人工粘性 収束するparticle pairに対してshock dissipationを与えます。 - `α`: linear viscosity - `β`: strong-shock penetrationを抑えるquadratic term shock捕獲に必要ですが、shear flowやdiskへ数値粘性を与えます。 ### 6.5 Ideal-gas EOS ```text P = (γ-1) ρ u c_s² = γP/ρ ``` `u`はspecific internal energyです。 ### 6.6 Tait EOS 液体demoではweakly compressible Tait EOSを使います。 ```text P = ρ0 c0²/γ_liq [(ρ/ρ0)^γ_liq - 1] ``` free surface付近の負圧は0へclipします。これはtensile instabilityを避ける教育用処理です。そのため圧力0の粒子が存在します。 ### 6.7 CFL timestep kernel lengthとsignal speedからglobal dtを選びます。強いshock、液体の高いsound speed、自己重力加速度によりdtは小さくなります。 ### 6.8 Smoothing length - fixed `h` - global adaptive `h` adaptive modeは全粒子共通のglobal `h`を平均neighbor数へ合わせます。particle-wise adaptive `h_i`やgrad-h correctionではありません。 ### 6.9 Self-gravity open boundaryではBarnes–Hutを用います。self-gravitating cloudとdiskは、SPH pressure forceとgravityを同時に時間発展させます。 ### 6.10 External potentials - uniform downward gravity: RT / dam break - central point mass: Keplerian disk / compact object - toy exponential cooling ## 7. ガス雲診断 ### Free-fall time 初期平均密度から ```text t_ff = sqrt(3π / 32Gρ_mean) ``` を計算し、`t/t_ff`を表示します。 ### Half-mass radius 全質量の50%を含む半径`R50`を計算します。 - `R50/R50,0 < 1`: 大域的収縮 - maximum densityだけ上昇: 局所集中またはnoiseの可能性 ### Radial velocity ```text v_r = r·v / |r| ``` 負は流入、正は流出です。 ## 8. Gas disk 中心点質量によるKeplerian circular speedを基準にします。 ```text v_phi ≈ sqrt(GM/r) Ω ≈ r^(-3/2) ``` 初期方位sectorをpassive labelとして保持し、差動回転を可視化します。 自己重力diskの`Q target`は初期sound speed設定用の目安です。局所surface densityとepicyclic frequencyから厳密に再計算したToomre Q mapではありません。 ## 9. Supernova remnant ### 9.1 Compact remnant prescription - illustrative non-monotonic islands - rapid-fallback toy 入力massとexplosion energyからneutron star / black hole、fallback、ejecta massを返します。これはhydrodynamic core-collapse solutionではありません。 ### 9.2 SPH ejecta 中心領域にradial velocityとhigh pressureを与え、ambient mediumへ膨張させます。表示するshock radiusは高圧particleのpercentileから得るproxyで、厳密shock finderではありません。 ## 10. 1D HLLC-ALE moving mesh ### 10.1 Conserved extensive quantities 各cellが保持するのは ```text M_i = ρ_i Δx_i P_i = ρ_i u_i Δx_i E_i = [P/(γ-1) + ρu²/2] Δx_i ``` です。 ### 10.2 Moving face flux face velocityを`w`とすると、moving faceを横切るfluxは ```text F_ALE = F(U) - wU ``` です。 ### 10.3 HLLC solver left / right acoustic waveに加え、contact waveを保持します。uniform pressure・velocityを持つdensity contactでは、meshがcontact speedへ追随すると相対mass fluxがほぼ0になり、固定meshよりadvection diffusionが小さくなります。 ### 10.4 Mesh motion interior face velocityは隣接cell速度の平均とregularization correctionから作ります。periodic contact testではdomain両端も同じbulk speedで移動し、表示時はdomain centerを差し引きます。 ### 10.5 Limitations - 1D - first order - fixed connectivity - no slope reconstruction - no limiter - no 2D/3D Voronoi tessellation - no topology change - no gravity coupling ## 11. Resource model ### Memory 推定には次を含みます。 - authoritative Float64 state - solver workspace - Worker copies - rendering frame - renderer overhead - diagnostic history - safety margin ### CPU / runtime solver complexity、particle / cell / mesh count、Worker数、推定step数から保守的に予測します。開始後の実測batchがより信頼できる値です。 ### Risk - Memory: RAM / JS heap fraction - CPU: Worker oversubscriptionとpredicted step time - Runtime: predicted total time - Overall: 最大risk ## 12. 数値収束の最低条件 1. dtを半分にする。 2. particle / cell数を増やす。 3. softening / smoothing / meshを変更する。 4. solver toleranceを厳しくする。 5. 主要診断量が一定値へ近づくか確認する。 6. seed依存性を確認する。 一つの画像だけではconvergenceを主張できません。