この文書は、現行 Fortran 実装の Coulomb FMM コア
bem_coulomb_fmm_core module page
と、その実装を分割した関連ファイル群の仕様とアルゴリズムをまとめたものです。
src/physics/field_solver/fmm/api/src/physics/field_solver/fmm/internal/common/src/physics/field_solver/fmm/internal/tree/src/physics/field_solver/fmm/internal/runtime/src/physics/field_solver/fmm/internal/periodic/対象は simulator 非依存の内部 API で、mesh_type や sim_config を直接 use しません。
BEACH 側では field solver adapter がこのコアを呼び出します。
この FMM コアの目的は、固定された source 点群 src_pos(3,n) と可変電荷 src_q(n) に対して、
多数の評価点で Coulomb 電場を高速に返すことです。
現行実装の設計目標は次の通りです。
free と periodic2 のみを対象にするコアが提供する主な手続きは次の 4 つです。
call build_plan(plan, src_pos, options)
call update_state(plan, state, src_q)
call eval_points(plan, state, target_pos, e)
call eval_point(plan, state, r, e)
入力・出力の意味は次の通りです。
src_pos(3,n):
source 点の座標。build_plan 後は固定とみなす。src_q(n):
source 点の電荷。update_state ごとに更新できる。target_pos(3,m) または r(3):
評価点。e(3,m) または e(3):
電場ベクトル。注意点:
k_coulomb を掛けていません。BEACH の adapter 側で最後に掛けます。build_plan は幾何依存処理、update_state は電荷依存処理です。eval_point(s) は plan と state が ready な前提です。BEACH の field solver adapter は、メッシュ要素重心を src_pos としてこのコアへ渡します。
build_plan の直後に update_state を実行します。plan を再利用し、src_q 更新として update_state だけを呼びます。build_plan と legacy tree metadata の同期をやり直すのは、plan 未構築時・source 数変更時・要素数 0 件で plan/state を破棄したときだけです。fmm_options_type主な内部オプション:
theta: well-separated 判定用パラメータleaf_max: source octree の葉に許す最大 source 数order: Cartesian 展開次数softening: softened Coulomb kernel の epsilonuse_periodic2: 2 周期軸モードの有効化periodic_axes(2), periodic_len(2): 周期軸と周期長periodic_image_layers: 近傍画像和の層数 Nperiodic_far_correction: core が受ける値は auto, none, m2l_root_oracle。periodic2 有効時の auto は m2l_root_oracle に正規化され、none は遠方補正なしとして残るperiodic_ewald_alpha, periodic_ewald_layers: m2l_root_oracle の build-time Ewald fit で使う分解パラメータと打切り深さtarget_box_min/max: dual-target tree を作るときの boxBEACH の adapter は現状 order = 4 を使いますが、コア自体は可変次数を受けられます。
periodic2 の auto は m2l_root_oracle に正規化されます。none は遠方補正を無効化します。
fmm_plan_type幾何にだけ依存する不変データです。
alpha, deriv_alphasource_leaf_nodesleaf_nodesnear_start/near_nodesfar_start/far_nodesm2l_target_nodes/m2l_source_nodesm2l_derivsource_p2m_basisfmm_state_type電荷に依存して毎回更新されるデータです。
src_q(n)multipole(ncoef, nnode)local(ncoef, n_target_nodes)multipole_active(nnode)local_active(n_target_nodes)multipole は source tree ノードごとの多重極係数、local は target tree ノードごとの局所展開係数です。
*_active は zero-node を早く飛ばすための 0/1 フラグです。
現行コアは softening 付き Coulomb kernel を使います。
近傍 direct 和でも遠方展開でも、同じ $G_\epsilon$ を使います。
多重指数 $\alpha = (\alpha_x, \alpha_y, \alpha_z)$ を使います。
node center を $c$ とすると、葉ノードの multipole 係数は
で定義します。
子ノード中心 $c_{\mathrm{child}}$ の係数を親中心 $c_{\mathrm{parent}}$ に平行移動して集約します。 $\mathbf{d} = c_{\mathrm{child}} - c_{\mathrm{parent}}$ とすると
現行実装では $\beta - \alpha$ に対応する index と
$\mathbf{d}^{\beta-\alpha} / (\beta-\alpha)!$ を build_plan 時に前計算します。
source node 中心 $c_s$、target node 中心 $c_t$ に対して $R = c_t - c_s$ とします。
局所展開係数は
で更新します。
ここで $D^\gamma$ は multi-index 微分です。
現行実装は $D^{\alpha+\beta} G_\epsilon(R)$ を m2l_deriv(:, pair) として pair ごとに前計算します。
親中心 $c_{\mathrm{parent}}$ の局所展開を子中心 $c_{\mathrm{child}}$ へ平行移動します。 $\mathbf{d} = c_{\mathrm{child}} - c_{\mathrm{parent}}$ とすると
これも build_plan 時に shift monomial を前計算します。
評価点 $x$ が属する target leaf の中心を $c_{\mathrm{leaf}}$ とし、 $\mathbf{dr} = x - c_{\mathrm{leaf}}$ とすると
で電場を評価します。 ここで $e_k$ は軸 $k$ の単位 multi-index です。
build_plan のアルゴリズムbuild_plan は幾何依存処理だけを行います。
source 座標の bounding box を再帰的に 8 分割して octree を作ります。 停止条件は次のどちらかです。
<= leaf_maxtarget 側は 2 通りあります。
target_box が無効:
source tree の葉をそのまま target leaf として使うtarget_box が有効:
box 全体を覆う別 target tree を作るperiodic2 では target point を box 内に wrap してから target leaf を探します。
各 target leaf について source tree を再帰走査し、 near node と far node を作ります。
well-separated 判定は
です。
free と periodic2periodic2 では $\mathbf{d}$ に minimum-image 補正を入れます。
その後、dual-tree 再帰で M2L pair cache を作り、 target node ごとの索引配列も準備します。
build_plan の最後で、refresh ごとに変わらない量を前計算します。
source_parent_ofparent_ofsource_p2m_basism2m_term_count, m2m_alpha_list, m2m_delta_listl2l_term_count, l2l_gamma_list, l2l_delta_listsource_shift_monomialtarget_shift_monomialshift_axis1, shift_axis2periodic_ewaldperiodic_root_operatorm2l_derivこの前計算により update_state は charge-dependent な加算だけに近づきます。
build_plan(src_pos, options):
initialize_basis_tables(order)
build_source_tree(src_pos)
precompute_source_p2m_basis()
build_target_topology(target_box)
build_interactions()
precompute_translation_operators()
precompute_periodic2_ewald_data()
precompute_periodic_root_operator()
precompute_m2l_derivatives()
update_state のアルゴリズムupdate_state は legacy 実装の refresh に相当する処理です。
source 座標は変わらず、src_q だけが変わる前提です。
update_state(plan, state, src_q):
ensure_state_capacity()
copy src_q
clear active flags
clear multipole/local only when the tree has no source leaves or no M2L pairs
P2M on source leaves
M2M bottom-up
M2L on cached pairs
L2L top-down
mark state ready
現行実装では、次の箇所に OpenMP を入れています。
update_state 全体を 1 つの parallel region で囲み、その内側で src_q コピーと active flag 初期化P2M: source leaf ごとのループM2M: 同一 depth の node ループM2L: target node ごとのループL2L: 同一 depth の node ループbuild_plan 時の translation / M2L 微分前計算各ループは「1 node 1 thread」になりやすいように書いてあり、 共有配列への更新は node 単位で独立させています。
update_state では次の無駄を避けています。
P2M の monomial 基底を source ごとに build 時に前計算するM2M/L2L の有効な (alpha, delta) 項だけを圧縮して持ち、無効項を走査しないM2L では source node の active flag を見て zero-node を pair 単位で早期 skip するM2L で target node 列へ細かく何度も書かず、thread-local な local_acc にためてから戻すP2M で target leaf ではなく source leaf 専用 index を使うeval_point(s) のアルゴリズム評価時の処理は次の通りです。
eval_point(r):
if plan is not built or state is not ready:
return zero vector
if periodic2:
wrap r into target box
leaf = locate_target_leaf(r)
if leaf not found or leaf is not mapped to a leaf slot:
use direct sum over all sources
if periodic2 and far correction is trunc:
add truncated far-image correction
if periodic2 and far correction is oracle:
add exact periodic Ewald correction
return
evaluate local expansion at leaf center
add near direct interactions
root local already carries periodic root correction when enabled
periodic2 では評価点を target box 内へ wrap してから探索するnear list に入った source index については direct 和を取ります。
periodic2 では [-N, N] x [-N, N] の画像シフトを陽に回します。
fallback でも同じ direct kernel を使いますが、periodic2 で m2l_root_oracle が有効なときは oracle 補正を別途加算します。
dual-target tree を使う場合、評価点が target box の外に出ることがあります。
そのときは target leaf を持たないので、全 source に対する direct 和へ fallback します。
m2l_root_oracle では build-time Ewald fit の teacher と同じ exact periodic correction を direct fallback へ足します。
m2l_root_oracle の root 補正は update_state 側で state%local(:, root) に注入されます。
したがって通常の leaf 評価では、eval_point(s) は root 補正を再計算せず、state に載っている local 展開をそのまま使います。
periodic2 と遠方補正periodic2periodic2 は「ちょうど 2 軸だけ周期、残り 1 軸は開放」です。
近傍画像和は
の有限画像を陽に足します。
M2L でも同じ画像シフト集合を使い、各 pair の derivative を画像和で前計算します。
periodic2 Ewald（Ewald2P）補正bem_coulomb_fmm_periodic_ewald.f90 は、2 周期・1 開放の Coulomb field に対する Ewald 形の補正を実装します。
ここでいう exact は「コードが実際に評価する有限和」を指します。理論上の無限和そのものではなく、field_periodic_image_layers = N と field_periodic_ewald_layers = L で real-space / reciprocal-space の打切り深さを決める build-time oracle です。
周期軸を a_1, a_2、開放軸を f とします。
周期長、セル面積、画像集合、逆格子集合を次のように置きます。
画像シフトと逆格子ベクトルは
と書けます。ソース位置を 、評価点を とし、
を導入します。以下では
field_periodic_ewald_alpha、
softening とします。
add_screened_point_charge が実装している screened Coulomb field は
です。これはポテンシャル
の勾配に一致します。
add_softened_point_charge が使う direct kernel は
で、通常の runtime direct path と同じ softening を使います。
実装上の real-space 補正は
です。実装では r2 <= tiny(1.0d0) の項はスキップするため、self interaction は入らない扱いです。add_periodic2_exact_ewald_correction_single_source に direct fallback の を足すと、inner image の softened 部分が打ち消され、outer shell 側は screened 形に置き換わります。
add_exact_periodic2_reciprocal_space_correction が使う逆空間項は、 に対して
を定義すると
です。コードでは term_p, term_m, pair_sum に対応します。
この式は、逆格子の k=0 を除いた高周波成分に対応します。
k=0 項add_exact_periodic2_k0_correction が実装しているゼロモード補正は
です。single-source の oracle では k=0 の電場寄与としてこの形を保持します。
以上をまとめると、add_periodic2_exact_ewald_correction_single_source が 1 粒子分に加える補正は
です。add_periodic2_exact_ewald_correction_all_sources はまずこれを全ソースに対して総和します。
charged_walls total-charge 補正非中性 slab の charged_walls closure に合わせて、add_periodic2_exact_ewald_correction_all_sources は全ソース和のあとに total-charge 補正
を加えます。ここで A = L_1 L_2 は周期セル面積、Q_tot = \sum_j q_j、z_low/high は target_box_min/max の非周期軸境界です。
この項は 2 枚の補償壁の場に対応し、slab 内では厳密に打ち消し合うため、target box 内で build する root oracle や通常の粒子前進には影響しません。影響するのは target box 外へ落ちた direct fallback 評価だけです。
field_periodic_ewald_alpha が <= 0 の場合、resolve_periodic2_ewald_alpha は
を自動選択します。min(L_1,L_2)\le 0 なら alpha = 0 として oracle を無効化します。
また内部では kmax = max(1, field_periodic_ewald_layers) として逆空間の有限和を作ります。
実際の runtime direct fallback は
です。m2l_root_oracle の build-time fit では check points が target box 内にあるため \mathbf E_{\mathrm{walls}} = 0 となり、teacher には single-source の \mathbf E_{\mathrm{corr}} だけが使われます。periodic_root_operator 側では定数 potential mode を使わないため、monopole column は 0 に固定されます。
m2l_root_oraclem2l_root_oracle は、この Ewald2P 補正を teacher にして root multipole から root local への演算子を proxy/check 点で fit するモードです。
periodic_image_layers = N: runtime で explicit に残す近傍画像殻periodic_ewald_layers = L: build-time oracle の real-space outer shell N < max(|i|,|j|) <= N+L と reciprocal cutoff |m|, |n| <= Lperiodic_ewald_alpha = alpha: Ewald 分解パラメータ。<= 0 なら自動決定local(:, root) += T_root_oracle * multipole(:, root) を加えるだけなので eval path に Ewald 本体は入らない固定次数 $p$、bounded interaction list を仮定すると、実用上は次のように見てよいです。
build_plan: $O(n \log n)$ に近いupdate_state: $O(n)$ に近いeval_point: $O(\log n + n_{\mathrm{near}} \, n_{\mathrm{img}}^2)$ に近いeval_points: 上記の点評価を target ごとに並列実行厳密な定数因子は次に強く依存します。
orderthetaleaf_maxperiodic_image_layersperiodic_ewald_layersこの FMM コアは汎用 kernel FMM ではありません。
order = 4build_plan 後に不変とみなすfree と periodic2periodic2 は正確に 2 周期軸が必要auto（既定）, none, m2l_root_oracle（periodic2 の default auto は oracle に正規化, none は補正なし）eval_point(s) の返り値には k_coulomb を含めない主な対応箇所:
src/physics/field_solver/fmm/api/bem_coulomb_fmm_core.f90,
src/physics/field_solver/fmm/api/bem_coulomb_fmm_core_build.f90,
src/physics/field_solver/fmm/api/bem_coulomb_fmm_core_state.f90,
src/physics/field_solver/fmm/api/bem_coulomb_fmm_core_eval.f90fmm_options_type, fmm_plan_type, fmm_state_type
（src/physics/field_solver/fmm/internal/common/bem_coulomb_fmm_types.f90）build_plan
（src/physics/field_solver/fmm/internal/tree/bem_coulomb_fmm_plan_ops.f90）update_state, p2m_leaf_moments, m2m_upward_pass, m2l_accumulate, l2l_downward_pass
（src/physics/field_solver/fmm/internal/runtime/bem_coulomb_fmm_state_ops.f90）eval_point, eval_points
（src/physics/field_solver/fmm/internal/runtime/bem_coulomb_fmm_eval_ops.f90）has_valid_target_box, use_periodic2_m2l_root_oracle,
use_periodic2_root_operator, build_periodic_shift_values, add_point_charge_images_field,
wrap_periodic2_point, apply_periodic2_minimum_image, distance_to_source_bbox,
distance_to_source_bbox_periodic
（src/physics/field_solver/fmm/internal/periodic/bem_coulomb_fmm_periodic.f90）resolve_periodic2_ewald_alpha, precompute_periodic2_ewald_data,
add_periodic2_exact_ewald_correction_single_source, add_periodic2_exact_ewald_correction_all_sources
（src/physics/field_solver/fmm/internal/periodic/bem_coulomb_fmm_periodic_ewald.f90）precompute_periodic_root_operator
（src/physics/field_solver/fmm/internal/periodic/bem_coulomb_fmm_periodic_root_ops.f90）src/physics/field_solver/bem_field_solver_config.f90,
src/physics/field_solver/bem_field_solver_tree.f90,
src/physics/field_solver/bem_field_solver_eval.f90設計上の責務分担は次の通りです。
mesh_type から src_pos を作る、q_elem を src_q へ流す、k_coulomb を最後に掛ける