BEACH アルゴリズム概要

この文書は、BEACH の現行 Fortran 実装に基づいて、数値アルゴリズムと実行順序を説明します。 BEACH は格子 PIC ではなく、三角形境界要素上の電荷を source とする BEM 風の Coulomb 場評価と、テスト粒子追跡を batch 単位で結合した表面帯電シミュレータです。

実装へのリンクは、現在のファイルと主な symbol を示します。行番号が後続変更でずれても、リンク先ファイルと symbol 名を優先してください。

アルゴリズム文書の構成

文書	内容
アルゴリズム概要	BEACH の計算モデル、初期化、batch ループ
場ソルバーと境界条件	Coulomb 場、direct/treecode/FMM 切替、periodic2 場境界
粒子追跡と表面電荷蓄積	粒子生成、Boris pusher、衝突、電荷堆積、統計・再開
Coulomb FMM コア詳細	FMM core API、tree 構築、M2L、periodic2 Ewald/oracle
`batch_duration` の安定性	batch 時間幅、定常値、線形安定性、Monte Carlo ノイズ

1. BEACH の計算モデル

Source: bem_simulator, bem_simulator_loop, bem_field_solver, bem_injection

BEACH の主状態は、三角形メッシュ要素の電荷 q_elem(i) と、batch ごとに生成される粒子群です。メッシュ要素は固定形状で、粒子は表面に衝突した時点で吸収されます。吸収された粒子の電荷は命中要素に蓄積し、次 batch の場計算へ反映されます。

1.1 状態変数

主要データ型は bem_types に定義されています。

型	主な内容	役割
`sim_config`	`dt`, `batch_count`, `max_step`, `field_solver`, `field_bc_mode`, `box_min/max`, `bc_low/high`	時間発展・場計算・境界条件
`mesh_type`	`v0/v1/v2`, `centers`, `normals`, `bb_min/max`, `q_elem`, `elem_surface_model`	三角形境界要素と電荷
`particles_soa`	`x`, `v`, `q`, `m`, `w`, `species_id`, `alive`	batch 内粒子の SoA 表現
`injection_state`	`macro_residual(:)`	`reservoir_face` の端数粒子数を batch 間・再開間で保持
`sim_stats`	`processed_particles`, `absorbed`, `escaped`, `survived_max_step`, `batches`, `last_rel_change`	実行統計

1.2 物理近似

各三角形要素は、重心 c_i にある点電荷 q_i として場に寄与します。要素面積上の連続電荷分布を積分する厳密 BEM ではなく、現行実装は重心点電荷近似です。

softening 付き Coulomb kernel:

G_\epsilon(\mathbf{r}) = \frac{1}{\sqrt{\lVert\mathbf{r}\rVert^2 + \epsilon^2}}

電位:

\phi(\mathbf{x}) = k_\mathrm{c} \sum_j q_j G_\epsilon(\mathbf{x} - \mathbf{c}_j)

電場:

\mathbf{E}(\mathbf{x}) = k_\mathrm{c} \sum_j q_j \frac{\mathbf{x} - \mathbf{c}_j} {\left(\lVert\mathbf{x} - \mathbf{c}_j\rVert^2 + \epsilon^2\right)^{3/2}}

ここで k_coulomb は bem_constants の Coulomb 定数です。粒子が実際に受ける電場は、境界要素電荷による場に一様外部電場 sim.e0 を加えたものです。

1.3 実行単位

BEACH は sim.batch_count まで batch を進めます。各 batch は次の意味を持ちます。

現在の表面電荷に基づき、粒子種ごとの粒子群を生成する。
現在の q_elem で field solver を更新する。
各粒子を最大 sim.max_step ステップまで追跡する。
衝突した粒子の電荷を要素別差分 dq に集計する。
dq を q_elem に加算し、必要なら表面モデルの緩和をかける。
統計・履歴を更新する。

sim.tol_rel は出力・監視用の値であり、現行 Fortran 実装では早期停止条件ではありません。

2. 実行エントリーポイントと初期化

Source: app/main.f90, bem_app_config_runtime, bem_mesh, bem_restart

2.1 main program の順序

app/main.f90 は CLI エントリーポイントです。概略順序は次の通りです。

順序	処理	主な実装
1	`--version` など設定読込前に完結する CLI option を処理	`handle_early_cli`
2	MPI と performance profile を初期化	`mpi_initialize`, `perf_configure_from_env`
3	設定読込、メッシュ構築、restart 読込または乱数 seed 初期化	`load_or_init_run_state`
4	履歴 CSV を open	`open_history_writer`, `open_potential_history_writer`
5	batch simulation を実行	`run_absorption_insulator`
6	summary と最終 CSV を出力	`print_run_summary`, `write_result_files`
7	RNG state と macro residual を checkpoint として出力	`write_rng_state_file`, `write_macro_residuals_file`
8	performance profile を出力し MPI を終了	`perf_write_outputs`, `mpi_shutdown`

設定ファイルは、明示引数があればその path、なければカレントディレクトリの beach.toml を使います。設定ファイルがない場合は default_app_config の既定値で走ります。

2.2 メッシュ構築

メッシュは mesh.mode に応じて template または OBJ から作られます。template の場合は build_template_mesh が次を行います。

mesh.templates を順に走査する。
enabled=false の template は飛ばす。
kind に応じて make_plane, make_box, make_cylinder, make_sphere などへ dispatch する。
template ごとに mesh_id を割り当てる。
surface_model と epsilon_r を要素配列へ展開する。
全 template の三角形配列を連結し、init_mesh へ渡す。

init_mesh は頂点配列から次を前計算します。

要素重心 centers(:, i)
要素法線 normals(:, i)
要素 AABB bb_min/max(:, i)
代表長 h_elem(i) = sqrt(area_i)
初期電荷 q_elem(i)
collision grid

2.3 periodic2 collision mesh

sim.field_bc_mode="periodic2" のとき、prepare_periodic2_collision_mesh は周期 2 軸に沿って各三角形を canonical な位置へ平行移動します。これは場計算の周期画像和とは別で、衝突判定用の primitive cell メッシュを安定化するための処理です。要素 index は base element のまま維持されるため、periodic image に命中しても電荷は base element に集約されます。

2.4 restart

output.resume=true のとき、load_restart_checkpoint は以下を読みます。

summary.txt: 完了済み batch と統計
charges.csv: 要素電荷
rng_state.txt または rng_state_rankNNNNN.txt: 乱数状態
macro_residuals.csv または rank 別 residual: reservoir の端数粒子数

sim.batch_count は累積到達 batch 数です。checkpoint が batches=100 で batch_count=150 なら、実行するのは 50 batch だけです。

3. batch ループ

Source: run_absorption_insulator, prepare_batch_state, process_particle_batch, commit_batch_charge

3.1 loop skeleton

run_absorption_insulator は、初期統計 initial_stats を受け取った場合はその batches から再開します。

final_batch_idx = sim.batch_count
batch_count_this_run = final_batch_idx - stats.batches

field_solver.init(mesh, sim)

for local_batch_idx = 1..batch_count_this_run:
    prepare_batch_state(...)
    field_solver.refresh(mesh)
    process_particle_batch(...)
    commit_batch_charge(...)
    count_batch_outcomes(...)
    MPI allreduce statistics
    accumulate_batch_stats(...)
    write charge/potential history when stride matches

3.2 batch state

prepare_batch_state は次を準備します。

今 batch の番号 batch_idx = stats%batches + 1
init_particle_batch_from_config による粒子群
dq_thread(nelem, nth): OpenMP thread ごとの電荷差分
photo_emission_dq(nelem): photo_raycast の放出元逆符号電荷
escaped_boundary_flag(:) と absorbed_flag(:)

dq_thread を thread 別に分けることで、衝突時の要素電荷加算を atomic なしで集計できます。

3.3 particle processing

process_particle_batch は粒子ごとに最大 sim.max_step 回の時間発展を行います。

順序	処理
1	現在位置 `x0` と速度 `v0` を読む
2	`field_solver%eval_e(mesh, x0, e)` で境界要素電荷による電場を評価
3	一様外部電場 `sim.e0` を加える
4	`boris_push` で候補位置 `x1` と候補速度 `v1` を計算
5	`find_first_hit(mesh, x0, x1, hit, sim=sim)` で線分衝突を調べる
6	hit があれば `q * w` を命中要素の `dq_thread(elem, tid)` へ加算し、粒子を吸収終了
7	hit がなければ `apply_box_boundary` で open / reflect / periodic を適用
8	粒子が生存していれば `x` と `v` を更新して次 step へ進む

BEACH_WARN_LONG_PARTICLE_STEPS を正整数で設定すると、長く生き残る粒子の診断出力を一定 step ごとに出します。

3.4 charge commit

commit_batch_charge は thread 別差分と photo emission 差分を統合します。

\Delta q_i = \sum_{\mathrm{thread}} \Delta q_{i,\mathrm{thread}} {}+ \Delta q_{i,\mathrm{photo}}

MPI 実行時は mpi_allreduce_sum_real_dp_array により rank 間で dq を和にします。その後:

q_i^{new} = q_i^{old} + \Delta q_i

を適用し、表面モデル緩和を実行します。最後に実際に変化した電荷量で監視値を計算します。

\operatorname{rel\_change} = \frac{\lVert\mathbf{q}^{\mathrm{new}} - \mathbf{q}^{\mathrm{old}}\rVert_2} {\max\left(\lVert\mathbf{q}^{\mathrm{new}}\rVert_2, q_\mathrm{floor}\right)}

この値は stats%last_rel_change と履歴出力に使われます。