🌊
【Unity/Compute Shader】内蔵GPUで動く！領域追従APIC-SWE-FFTハイブリッド流体　後編

2026/03/16に公開
 4. 既存の流体シミュレーション手法の比較と限界海全体を表現するにあたり、単一のインタラクティブな計算手法で構築することは最新のGPUをもってしても難しいと思います。今回は限られたリソース(内蔵GPU)での動作を目標としているため、プレイヤーからの距離に応じて計算手法を切り替えるLOD (Level of Detail)の導入が必須 となります。
手法を一つずつ解説すると相当長くなってしまうため、ここでは、今回比較検討のために実装した代表的な流体シミュレーション手法について、2次元断面を載せつつ要点のみ解説します。
流体は一般にナビエ・ストークス方程式に従いますが、これを数値的に解くCFD(数値流体力学)のアプローチは、空間を固定して観察するオイラー的な見方(格子法)と、流体を粒子の集まりとして追いかけるラグランジュ的見方(粒子法)、そしてそのハイブリッド手法に大別されます。


格子法のイメージ
粒子法のイメージ






 各手法の特徴と評価
 MAC (Marker and Cell)
アプローチ: 格子法
概要: 空間を固定のグリッドに分割し、そこに圧力や速度を保持します。速度を面上で管理し、圧力をセル中心で管理する「スタッガード格子」を用いることで、圧力計算が安定します。
性質: 水面の表現は比較的綺麗に出ますが、移流(流体が動く処理)をグリッド上で行う性質上、数値的な拡散が起きやすく、激しい水しぶきなどの表現には不向きです。


スタッガード格子(MAC格子)

 SPH (Smoothed Particle Hydrodynamics)
アプローチ: 粒子法
概要: 流体をすべて粒子として扱い、周囲の粒子との距離に応じた重み付け（カーネル関数）で物理量を計算します。
性質: 質量保存が確実で、水しぶきや自由表面の表現が容易です。しかし、圧力を安定させるのが難しく、圧縮性のあるプルプルとした挙動になりがちです。また、多くの粒子を用いるため、計算量が膨大になります。

 PIC (Particle-in-Cell)
アプローチ: ハイブリッド (格子法 ＋ 粒子法)
概要: 格子法と粒子法のいいとこ取りをした手法で、移流を粒子で計算(Advection) → 格子に当てはめ(P2G) → 計算の難しい圧力の反復計算を格子で行う(Pressure Solve) → 粒子に速度を書き戻す(G2P)という流れで行います。
性質: グリッドと粒子の間で速度データを補間してやり取りする際、エネルギーが大きく散逸してしまいます。結果として、スライムのようなもっさりした動きになってしまいます。


PIC法の流れ

 FLIP (Fluid-Implicit-Particle)
アプローチ: ハイブリッド (格子法 ＋ 粒子法)
概要: PICの「粘り気が出る」という弱点を克服するため、速度そのものではなく「速度の変化量」を粒子に書き戻す手法です。安定のため、数パーセントPIC法を混ぜるのが一般的です。
性質: PICよりも圧倒的にダイナミックで水らしい動きになりますが、今度は逆にノイズが乗りやすく、表面がボコボコと暴れやすくなるという別のトレードオフが発生します。計算量はPICの1.1～1.3倍程度です。

 APIC (Affine Particle-in-Cell)
アプローチ: ハイブリッド (格子法 ＋ 粒子法)
概要: 粒子が単なる速度ベクトルだけでなく、局所的な速度の微分(アフィン行列)を保持する手法です。現在、HoudiniやBlenderなどのDCCツールにおいてFLIP法とともに主流となっているアルゴリズムです。
性質: 粒子が速度だけでなく局所速度勾配（affine velocity field）を保持し、渦や飛沫の形状などを美しく計算できる手法です。PICほどではないですが、FLIPと比べると多少粘性が強くなります。計算量はFLIPの1.3～1.8倍程度です。

 SWE (浅水方程式 / Shallow Water Equations)
アプローチ: 2次元近似法
概要: 「水深が水平方向の広がりに比べて十分に浅い」という仮定のもと、鉛直方向を平均化した深さ平均方程式です。(上の例は1次元)
性質: 3次元の計算を2次元の高さマップとして扱えるため、圧倒的に軽量です。波の伝播や障害物とのインタラクションもそれらしく表現できますが、波が巻き込んで砕けるような三次元的挙動や、飛沫を表現することは単体ではできません。


SWEの表現の限界

 FFT (高速フーリエ変換 / Fast Fourier Transform)
アプローチ: 統計的波形合成
概要: 厳密には流体ソルバではなく、海洋学の統計データに基づいた波の合成です。
性質: 広大な海域でもGPU上で非常に高速に計算できるため、背景の海として有用です。しかし、あくまで「見た目だけ」の波であるため、キャラクターが飛び込んだり、船が波をかき分けたりするインタラクションはできません。

 手法まとめ


FFT
SWE
MAC
SPH
PIC
FLIP
APIC


インタラクション
×
〇
〇
〇
〇
〇
〇

波の巻き込み表現
×
×
△
〇
〇
△
〇

安定性
◎
〇
◎
△
〇
△
〇

数値粘性の低さ
粘性なし
△
△
△
×
〇
〇

GPU計算効率
◎
〇
△
×
△
△
FLIPの約1.5倍

特徴
触れない
1次元低い
格子法
粒子法
格子＋粒子
PIC改良
FLIP改良

※数値粘性が低いほど数値拡散によるエネルギー損失が小さい

 5.手法の選定当初掲げた「見た目の自然さ ＞ 軽量さ ＞ 物理的な正しさ ＞ 実装の簡単さ」という選定基準と、Pythonでの事前プロトタイピングによる検証を経て、最終的にAPIC、SWE、FFTの3手法をLOD（Level of Detail）的に組み合わせる構成に決定しました。
それぞれの弱点を補い合う設計思想は以下の通りです。
APIC: 水しぶきや巻き込みといった立体的なインタラクション表現において最高峰のルックを持ちますが、計算負荷がリアルタイム用途では非常に高いという課題があります。そのため、プレイヤーの直近（最も激しく水が動く領域）かつ、必要な生存粒子のみに限定して稼働させます。
SWE: 負荷が低く、波紋の伝播を綺麗に表現できますが、立体的な波は作れません。そこで、APICの外側を取り囲むインタラクション領域として配置し、波が激しくなった瞬間だけ局所的にAPICへ処理を引き渡す(昇華させる)土台として運用します。
FFT: 負荷は非常に低いですが、干渉できません。そのため、プレイヤーから遠く離れた背景の海として配置します。

 各手法の構成と採用アルゴリズム各流体ソルバをUnityのCompute Shaderに落とし込むにあたり、リアルタイム動作の要となる計算手法には、描画速度と安定性を最優先したアルゴリズムを選定しています。

 APICの内部構成圧力計算：Jacobi-PCG法(前処理付き共役勾配法)

非圧縮性流体では圧力のポアソン方程式を解くステップが必要ですが、十分解に近づくまで何周も反復(iteration)する必要があり、通常、最も計算負荷が高くなるステップです。

これを解く一般的な手法に、Jacobi法、CG法、PCG法、MultiGrid-CG法などがあります。MultiGridは領域の大きい場合でも計算量が増えづらいのが利点ですが、狭い範囲では逆に重くなってしまいます。REIMUではローエンド環境で最速と思われる(Jacobi-)PCG法を採用しています。
メッシュ化：Marching Cubes

生存しているAPIC粒子から密度のスカラー場を生成し、リアルタイムにマーチングキューブでポリゴン化しています。

生成した頂点データはCPUに読み戻さず、DrawProceduralIndirect を用いてGPU完結でフルードメッシュや飛沫の描画を行うことで、転送ボトルネックを排除しています。さらに、少し離れた点の密度(空気混入度)を活用することでレイトレーシングなどの重い処理を避けつつ泡領域と水領域に分けています。

 SWEの内部構成計算スキーム：有限体積法(FVM), HLL近似リーマンソルバ

SWE(浅水方程式)のような双曲型偏微分方程式を解く際、単純な差分法では波がすぐに減衰してのっぺりしてしまったり、逆に数値振動で破綻(爆発)したりします。波の鋭さを保ちつつ安定させるため、流束の計算にHLL近似リーマンソルバを用いています。
高解像度化：MUSCLスキーム ＋ MinModリミッター

さらに波面のディティールをパキッと出すため、空間精度を上げるMUSCLスキームを導入しました。MinModリミッターで極値をカットすることで、どんなに激しくかき混ぜても爆発しない堅牢な水面を実現しています。
泡の移流： SWEの速度場を利用して、発生した泡(白波)のテクスチャ座標を移流シミュレーションさせています。
Dry/Wet境界の安定化： 極端な浅瀬での数値爆発(ゼロ除算)を防ぐため、マニングの粗度係数を用いた底面摩擦の適用や、地形の壁によるフラックスの遮断など、堅牢な境界処理を実装しています。

 FFTの内部構成スペクトル：Jonswapスペクトル

海洋工学などで用いられる、風が吹き荒れる海の波のスペクトルモデルです。
LODとモーフィング：

背景の海を3段階のカスケード(解像度の異なるメッシュ)で構成しています。単純に切り替えると境界のLODポップ(メッシュの急な切り替わり)が目立つため、頂点シェーダー側でカメラ距離に応じたモーフィングを行い、シームレスにブレンドしています。

 6. 異なる流体を繋ぐ遷移の実装(Coupling)REIMUの最大のポイントは、APIC(3D粒子)とSWE(2D高さマップ)という次元の異なる流体シミュレーションを、破綻なく、かつリアルタイムに双方向で遷移させている点にあります。

 通常のLODからの脱却当初、APICの領域の外側にSWEを配置し、「水平方向」で流体を繋ごうと試みました。間に緩衝地帯のようなものを設けて、そこでAPIC粒子の増減を行い、それをSWEに高さとして伝えるというものです。しかし、領域の境界でしか粒子の増減が発生しないため、APIC粒子が中央からどんどん流れ出てしまったり、緩衝地帯で新たに生まれた粒子が流れていかずに詰まってしまったりしました。

さらに、この方法ではAPIC領域が固定なため、静かな波であっても計算負荷が大きくなってしまいます。
この2点を一気に解決するために導入したのが、鉛直方向に領域を重ね、SWEをAPICの「床」として扱うというアプローチです。このSWE-APIC間の遷移は「波の激しさ」(後述)に応じて自動的に行われるようにしています。「床」と表現しましたが、常に2層あるわけではなく、SWE100%で計算負荷抑えたり、APIC100%で水が吹き飛ばされたあとの海底を見せたりすることも可能です。

これにより、SWEで表現しきれないときだけ、その場所だけ APICに切り替える、ということが可能になり、計算負荷が大幅に下がります。

 SWE → APIC（昇華：Sublimation）SWEが上方向に向かってAPIC粒子をスポーンさせます。これを「昇華」と呼ぶことにします。

 昇華判定3つの物理量から閾値を設定して判定しています。
フルード数 (Froude Number): 流速と浅水波速度の比。Fr = \frac{|\mathbf{v}|}{\sqrt{gh}}
高さのラプラシアン: 水面の凹凸・尖り具合。\nabla^2 h
水面の勾配 (Gradient): 波の急峻さ。\nabla h
これらの値が閾値を超えた箇所のSWEセルを、APICの発生源とみなします。

計算負荷を下げるため、ルート計算を避けてフルード数の2乗で判定するなどの最適化を行っています。

 鉛直方向への運動量変換SWEは2次元の計算手法なので、鉛直方向の速度をそもそも保持していません。そこで、非圧縮性流体の連続の式 \nabla \cdot \mathbf{u} = 0 に着目しました。水平方向の速度の変化（集まり具合）は、鉛直方向の速度（湧き上がり）に変換されるはずです。REIMUでは、SWEから取得した高さのラプラシアンのマイナス成分(水面が凹から凸へ急激に盛り上がる力)を、疑似的に鉛直方向の初速度 push_z として与えています。これにより、波がぶつかり合って「上に逃げ場を失った水」が、非常に自然な水しぶきとして空中に舞い上がります。
// 実際のCompute Shaderにおける初速度の計算
float c_wave = sqrt(9.81f * safe_h);
// ラプラシアンの負の成分（波の盛り上がり）を上方向の速度に変換
float push_z = clamp(max(0.0f, -laplacian_Eta) * push_z_multiplier, 0.0f, c_wave);
p.velocity.z = push_z + (rz * splash_velocity_multiplier) * c_wave; 

 APIC → SWE（凝縮：Condensation）空中に舞い上がったAPIC粒子が海面(SWEの表面)に落下してきた際、粒子を消滅させ、その質量と運動量をSWE側に書き戻します。これを「凝縮」と呼ぶことにします。

 凝縮判定SWEの下にAPIC粒子が潜ったら凝縮させます。

 運動量の加算落下してきた粒子をただSWEの水位に足すだけでは、波紋がきれいに広がりません。そこで、落下したAPIC粒子の鉛直方向の運動量を、SWEの水平方向の運動量
(hu, hv)に変換して放射状に拡散させています。
// 落下の衝撃（垂直方向の運動量）
float impact_energy = vertical_momentum_to_wave * p.mass * abs(p.velocity.z) * weight;

// 衝撃を放射状の波（水平方向の運動量）に変換
float delta_hu_wave = impact_energy * radial_dir.x * depth_factor;
float delta_hv_wave = impact_energy * radial_dir.y * depth_factor;

 保存則と安定化の工夫SWEとAPICを同一流体と捉えて変分することによって、厳密に保存則(質量・運動量・エネルギー)を守ろうとしたこともありましたが、それでは見た目が破綻してしまいます。
事前に「見た目の自然さ ＞ 軽量さ ＞ 物理的な正しさ ＞ 実装の簡単さ」という優先順位は決めていたので、厳密な実装は研究に譲り、REIMUはゲーム用と割り切って、
質量は保存(数値計算上の誤差は許容)
運動量は見た目を損なわない範囲で保存
エネルギーは安定化のため意図的に散逸(現実でも音や熱となる)
という方針を定めました。
さらに、安定化のために
粒子に寿命を持たせ、昇華してすぐのAPIC粒子は凝縮させない
1フレームにSWEからAPICに昇華できる高さに限界を設ける
凝縮時に周囲のSWEセルにも分配する
といった工夫も行っています。

 SWE ↔ FFT間の遷移当初、SWEの境界条件としてFFTを与える方針でしたが、波が尖っているFFTに比べ、SWEが滑らかすぎたので、SWEにFFTを高さとして加算する方式にしました。これにより、通常インタラクション不可なFFTの波を直接触れるという体験を提供しています。

 7.Compute Shaderのアーキテクチャとデータフロー実装はUnityのCompute Shaderで行いました。Compute Shaderとは、GPUを用いて並列計算を効率的に行うもので、"Shader"から連想される「影付け」だけでなく、汎用的な計算が可能です。
UnityではCompute Shader側でHLSL言語で記述したカーネル（関数）に対して、並列のスレッド数を指定し、C#のManager側からDispatchという命令を送ることで動作します。情報の書き出し・読み出しは、カーネルにセットしたBufferやTextureで行います。
GPUは同じ動きを全体でやる並列計算に向いており、1つの複雑な計算を行うCPUに比べて圧倒的な計算速度を誇ります。 一方で、1つのセルに対して複数のスレッドから情報を書き込む 「アトミック加算」がボトルネック になりやすく、このような処理をできるだけ減らすことが最適化に繋がります。
例えば圧力計算（PCG法）で必須となる巨大な配列の内積計算では、groupshared メモリを用いた並列リダクション(Parallel Reduction)を実装し、アトミック加算の回数を劇的に減らしています。

さらに、全粒子ではなく 「生きているAPIC粒子」に絞って計算を行うことで、一番重い処理である圧力の反復計算やボクセル化におけるアトミック加算の負荷を最小限に抑えています。

 ファイル構成とデータ・命令の流れ

yEdで作成

 ディスパッチ（Dispatch）のフロー地形の取得
マウス入力
スレッド数の最適化
FFT: スペクトル計算
FFT: バタフライ演算
APIC: 移流 (Advect)
APIC → SWE: 凝縮 (Condense)
SWE: Solve
SWE → APIC: 昇華 (Sublimate)
APIC: P2G (Particle to Grid/ 粒子を格子に当てはめ)
APIC: Jacobi前処理
APIC: 圧力ポアソン反復
APIC: 密度計算
APIC: G2P (Grid to Particle/ 格子の情報を粒子に書き戻し)
APIC: ボクセル化
APICボクセル, APIC飛沫, SWE・FFTの描画命令の送信

 8.性能
 検証環境OS: Windows11
CPU: Intel Core i5-8500
GPU: Intel Graphics 630
RAM: 16GB
Unityバージョン: Unity6 (6000.3.6f1 LTS)
描画解像度: 1920×1080

 設定デフォルトの設定は以下です。

*APIC: 64^3, 128^3でボクセル化

*SWE: 128^2

*FFT: 512^2, 3段階カスケード

 動作速度2017年モデルのiGPU (Intel Graphics 630)を用いたFHD解像度での厳しいテスト環境においても、波を立てていないときは約83ms (約12FPS)、マウスで激しくかき混ぜたときでも約106ms (約9FPS)で動作します。

下の表からもわかる通り、ボトルネックは「APICのボクセル化」と「透明オブジェクト（水面）の描画」に集中しています。流体シミュレーションの純粋な計算部分（移流や圧力計算など）だけで見れば約50msに収まっており、近年のディスクリートGPU環境であれば余裕をもって60FPS動作が狙えるアーキテクチャになっています。

 フレームタイム内訳ビルド後、RenderDocで測定。


カテゴリ
処理名
所要時間


Compute
SWEシミュレーション
1.3ms


APIC 圧力計算
5.2ms


ボクセル化
29.2ms


その他(移流,P2G,G2Pなど)
45.1ms

Render
不透明オブジェクト
4.6ms


透明オブジェクト
40.3ms


その他(SSAO等)
16.3ms

合計

106.3ms


 既存のソルバとの比較同じ環境で既存の手法と比較しました。


手法
最小処理時間(※1)
最大処理時間(※2)
備考


フルAPIC
約7,700ms(※3)
約8,900ms(※3)
FFT領域まで全てAPICで計算した場合の推定

APIC+FFT
約168ms
約186ms
昇華・凝縮のシステムがない場合の負荷

REIMU(本作)
約83ms
約106ms
APIC-SWE-FFTカップリング

SWE+FFT
約78ms
約91ms
APICの立体的な波を諦めた場合の負荷

フルFFT
約47ms
約47ms
インタラクションを諦めた場合の負荷

※1 そっと水の中に入った状態で1フレームの処理時間をそれぞれ数回測定し、最も時間が短かったもの。
※2 できるだけ激しくかき混ぜた状態で1フレームの処理時間をそれぞれ数回測定し、最も時間が長かったもの。
※3 これのみ推定値。APICで64^3を計算するのにかかる時間はおよそ
 \rm{64^3APIC}\approx\rm{(64^3APIC+512^2FFT)}-\rm{(512^2FFT)}APICで512\times 512 \times 64を計算する場合、グリッド数が64^3の64倍なので、計算量もおよそこの64倍。

 9.技術的課題と解決策(FAQ)開発中、行き詰まった点に対するフェイクや最適化での解決策をQ&A形式でまとめました。
Q. APICとSWEを水平方向で繋ぐと、境界に粒子が溜まってしまう

A. 領域を重ね、鉛直方向(上下)に繋ぐ構造へ変更

水平に並べて繋ごうとすると次元の違いから不自然な壁ができてしまったため、APICとSWEの領域を完全に重ね合わせ、水面から粒子が飛び出す「昇華」と、水面に落ちて波紋になる「凝縮」という上下のやり取りで解決しました。
Q. GPUの浮動小数点アトミック加算(InterlockedAdd)が重い

A. 固定小数点ハックによる最適化

InterlockedAdd は通常、整数にしか対応しておらず、浮動小数点でCASループを回すと著しく重くなります。そこで、加算前に整数型にキャストして加算し、読み出し時に元に戻す疑似的な固定小数点演算を挟むことで、精度を保ちつつ高速化しました。
Q. APICの計算負荷が高すぎてリアルタイムで動かない

A. 生存している粒子だけに計算を絞る

最大粒子数で常に計算を回すのではなく、毎フレーム「生きている粒子」だけのインデックスリストを構築し、そこに対してだけDispatchを行うようにして、最も重い圧力計算の負荷を劇的に下げました。
Q. リアルタイムのコースティクスだと、SWEとFFTの境界線が露骨にバレる

A. テクスチャを水面法線で歪ませるフェイク手法へ

厳密な光の屈折シミュレーションは諦め、事前に手描きしたコースティクステクスチャを、水面側の法線を使って歪ませるだけのシンプルなルックデヴに落とし込みました。
Q. 地面側にコースティクスを描画すると、陸の上(水がない場所)にはみ出る

A. 海面側のシェーダーから深度を参照して描画

地面のシェーダーに処理を書くと、波が激しいときにどうしてもコースティクスが見えてしまいました。そこで、「水面のシェーダー」側でコースティクスが水底に張り付いて見えるように描画することにしました。
Q. APICボクセルの表面で密度が一定になってしまい、泡の情報が取得できない

A. 少し離れた位置の密度を用いる

APICの密度については以下のような手法を検討しましたが、うまくいきませんでした。

・APICは全部泡にする→SWEとの切れ目が目立つ

・APIC粒子の下側にあるSWEの泡情報を参照する→不自然な切れ目が発生

・ボクセルの(絶対)高さを利用→同じ水滴なのに、高度が下がると急に水戻って変

・APICの密度を利用→ボクセル化されるほど密度がある場所では表面の密度がすべて一定になってしまう。

そこで、「ボクセル化or飛沫判定」用の密度とは別に、各セルに対して「周囲2セル離れた位置の密度を6点(正規化してから)サンプリングしてきた平均」を計算することで、対応しています。この際、水が薄い部分で地形の下を読み取ってしまっていたので、地形の下はすべて水判定で計算しています。
Q. APICの飛沫を球状のビルボードで描画すると違和感がでる

A. 形状をひし形にして、速度に応じて引き延ばし、ランダムな発光のみに

実際は飛沫は円または楕円が正しいはずなのですが、三角形や四角形で描画されているSWEやFFTと並ぶと球は悪目立ちします。そこで、形状をひし形に変更し、速度に応じて引き延ばすようにしました。さらに、SWEがちょうど太陽を反射したときにだけキラッと光る現象に着目し、APICの飛沫も「普段は透明にし、ランダムなタイミングで発光させる」ようにしたことで、どれだけカメラの近くに飛んできても違和感なく見えるようにしました。
Q. 真上から撮影する正投影手法だと、地形の高さが上手く取得できない

A. レイキャストによるハイトマップ取得へ変更

オープンワールドのように無限に広い海ではないため、初期化時にC#側で地形に向けてレイキャストを飛ばし、その結果をハイトマップテクスチャに焼き付けてCompute Shaderに渡す堅実な方法を取りました。

 10. おわりに
 AIを活用したモダンな開発フローAPICは2015年の論文で発表された比較的新しい方法で、ネット上の実装例がまだ少なく、AIの情報の精度がかなり悪いです。そのため、事前にPythonでテストを重ねてプロトタイプと仕様書を作成し、それをもとにAIをディレクションしてC#/ Compute Shaderに落とし込みました。これにより、Unityをインストールし、構想を開始してから1か月半ほどで完成まで持っていくことができました。
当初から、プレイヤーの目から近くて違和感のないカップリングが求められるAPIC-SWE間が鬼門だと思っていたので、Pythonで以下のような試作を繰り返し、設計が固まった段階で実際に作り始めました。



FLIPとSWEのカップリングテスト

合計質量(FLIP粒子数＋SWE面積×密度)が保存している
この後に行ったバグ修正・仕様追加の影響で完成版と仕様が異なる箇所もありますが、実際に使用した当初の仕様書(プロンプト)が以下です。
仕様書ハイブリッド流体シミュレーション (2D SWE + 3D APIP) コア仕様書

0. --- 前提 ---
本プロジェクトの大目標は「海全体を高精度で描画すること」である。
優先順位は常に、「視覚的な自然さ(違和感のなさ)」→「プログラムの軽量さ」→「物理的な正しさ」→「実装の簡単さ」とする。

1. --- 概要 ---
広大な海域全体を均一な高解像度でシミュレーションすることは不可能なため、低負荷な2次元SWE（浅水方程式）と高精細な3次元APIC（Affine Particle-In-Cell）を動的LODとしてシームレスに結合する。
基本構成として「SWEの動的な水面を境界（土台）とし、その上でAPICの層が動く」形をとる。砕波などの激しい現象が起きる領域では、水柱全体がAPIC粒子へと「昇華」し（床から水面まで丸々APIC化）、波が穏やかになった領域ではAPIC粒子がSWEへと「凝縮」することで、シームレスなマルチスケール流体を実現する。

2. --- 空間と質量の定義 ---
計算は論文に合わせてZ-upで行い、描画時にUnityに合わせてY-upに変換する。

・低次モード: 2次元SWE（有限体積法）
x-y平面上の2次元グリッド（セル幅 dx_swe, dy_swe。基本は正方形 dx_swe = dy_swe）。各セルが水深 h と水平流速 u, v を持つ。
セルの質量：M_swe = rho * h * dx_swe * dy_swe
保存変数：U = (h, hu, hv)
x方向フラックス：F(U) = (hu, h*u*u + 0.5*g*h*h, h*u*v)
y方向フラックス：G(U) = (hv, h*v*u, h*v*v + 0.5*g*h*h)

・高次モード: 3次元APIC
x-y-z空間を運動する粒子群。各粒子は位置(x, y, z)、速度(u, v, w)、質量 m_p、および3x3のアフィン行列 Cp を持つ。
粒子のベース質量 m_base は固定定数とし、生成時はこの質量を基準に粒子数を決定する。

・解像度比の制約
「dx_swe = M * dx_apic」とする。M は必ず2以上の整数。APICからSWEへのインデックス変換は「i_swe = floor(x_apic / dx_swe)」「j_swe = floor(y_apic / dy_swe)」で確定させる。

・保存の原則
全質量は常に一定。

3. --- 実装するプラットフォーム ---
・エンジン: Unity（URP）
・制御層: C# (マネージャー)
・演算層: HLSL (Compute Shader)
・アトミック演算の仕様: SWEへの書き込みにおいて、激しいスレッド競合（Spin-lock状態）によるGPUフリーズを防ぐため、直接SWEの現行ステートを加算更新することは禁止する。必ず中間ソースバッファ（Delta_H_Buffer, Delta_HU_Buffer等）を用意し、そこにCASループでアトミック加算を行った後、別カーネルでSWEステートに反映させる。

4. --- 実装の手順 ---
Step 0: 動的 dt とサブサイクリング回数の決定
SWEは高速な波を扱うため、APICの1ステップ(dt_apic)の間に N回 のサブサイクリング(dt_swe)を行う。
wave_speed_max = 領域内最大値( sqrt(u*u + v*v) + sqrt(g * max(h, h_min)) )
dt_swe_max = C_cfl * dx_swe / wave_speed_max （C_cfl はもっさり感を防ぐため限界の0.8程度に設定）
N = ceil(dt_apic / dt_swe_max)
dt_swe = dt_apic / N

Step 1: APIC移流・外力適用
全粒子に対して以下を行う。
・重力を鉛直速度 w に加算する。
・3x3アフィン行列 Cp を用いた移流により位置と速度を更新する。
・壁や地形との衝突処理を行う。

Step 2: 凝縮フェーズ（APIC → SWE）
目的は、不要粒子の削減と大域水面への滑らかな復帰である。単一セルへの集中投下による「無限の柱」や「鋭い棘」の発生を防ぐため、スプラッティング（周辺セルへの分散加算）を行う。
・対象条件: z <= h に存在する粒子は強制的に凝縮対象とする（SWE領域内部への粒子の潜り込み・取り残しをアルゴリズム的に許容しない）。h は粒子の水平位置に基づきバイリニア補間して取得する。ただし、生成直後の点滅を防ぐため、生成後一定時間（例: 0.1秒）経過していない粒子は除外する。
・運動量変換：
粒子の水平運動量：delta_hu_mass = m_p * u_p * weight、delta_hv_mass = m_p * v_p * weight
粒子の鉛直運動量を水平方向へ分配（波のチャタリング防止のためデッドゾーンと正規化勾配を使用）：
grad_norm_x = grad_h_x / (abs(grad_h_x) + epsilon)
grad_norm_y = grad_h_y / (abs(grad_h_y) + epsilon)
delta_hu_wave = gamma * m_p * abs(w_p) * grad_norm_x * weight_smooth
delta_hv_wave = gamma * m_p * abs(w_p) * grad_norm_y * weight_smooth
・SWE中間ソースバッファへの加算処理:
落下地点を中心とした 3x3 または 5x5 のセルに対し、ガウスカーネル等の重み(splat_weight)を用いて分散させる。
質量ソース: delta_h_i = (m_p * splat_weight) / (rho * dx_swe * dy_swe)
運動量ソースu: delta_hu_i = (delta_hu_mass + delta_hu_wave) * splat_weight / (rho * dx_swe * dy_swe)
運動量ソースv: delta_hv_i = (delta_hv_mass + delta_hv_wave) * splat_weight / (rho * dx_swe * dy_swe)
これらを Delta_H_Buffer 等へアトミック加算する。
・処理後、粒子を削除する。

Step 3: SWEサブサイクリング
APICの1ステップ中に、SWEを N回 更新する。
source_step = Intermediate_Buffer_Value / N

・Step 2で中間ソースバッファに蓄積された値を source_step として現行のSWEステートに加算する。
・HLLCによる2次元移流計算を行う。
・ドライステート問題（底到達時の暴れ）を防ぐため、h < h_min のセルは計算を完全にスキップし、速度を0にリセットする。流速計算時はゼロ除算を避け、 u = (sqrt(2) * h * (hu)) / sqrt(h^4 + max(h, epsilon)^4) などの非特異化を行う。
・TVDリミッターにより安定化する（波の鋭さを保つためSuperbee等の検討も可）。
・微小な振動を吸収するため、底面摩擦項（マニングの粗度係数など）を運動量にわずかに適用し、エネルギーを散逸させる。

Step 4: 昇華フェーズ（SWE → APIC)
目的は、砕波や激しい変形部の局所高解像度化である。砂山状態を防ぐため、表面を削るのではなく「激しい部分は床から水面まで丸々APICの層」にする水柱全体変換を行う。
・対象条件:
Fr平滑化
Fr = sqrt(u*u + v*v) / sqrt(g * max(h, h_min))
Fr_smoothは周囲セルを用いた2Dガウシアンフィルタ
Fr_filtered = lambda * Fr_prev + (1-lambda) * Fr_smooth
・生成量（水柱全体変換）:
Fr_filtered > Fr_threshold を満たした場合、そのセルの h を h_min まで一気に減算し、減少した体積 (h - h_min) 分をすべてAPIC粒子に変換する。
・粒子生成：
生成すべき総質量 = rho * (h - h_min) * dx_swe * dy_swe
生成粒子数 = floor(生成すべき総質量 / m_base)
各粒子の質量 = 生成すべき総質量 / 生成粒子数（または m_base 固定とし端数はSWEに残す）
x, y配置: セル内にランダムまたは均等分布
z配置： [h_min, h] 内に分布（床から水面までの柱状に配置）
u, v = SWEのu, v
w = dh_dt + 初期上向きノイズ（生成直後の即座な凝縮・点滅を防ぐためのスプラッシュ初速）
Cp = ゼロ行列で初期化（生成後周囲の流れに自然に馴染ませる）

Step 5: APICグリッド転写（P2G）と動的境界条件
・粒子の質量と運動量をAPIC背景3Dグリッドへ転写する。
・境界振動の平滑化フィルタ（床の暴れ防止）: SWEの細かい波紋がAPICを激しく叩き上げるのを防ぐため、水面高さ h と鉛直速度 dh_dt に対して強力なローパスフィルタ（空間的・時間的な平滑化）をかけたものをAPIC側の境界として使用する。
・動的剛体境界の滑らかな評価（階段状ノイズの排除）:
APICグリッドの圧力を計算する際、セルの状態を SOLID(0) / FLUID(1) の2値で判定せず、「セルが水に浸かっている割合（Fractional Volume）」 phi を計算する。
phi = clamp((h - (z_k - 0.5 * dz)) / dz, 0.0, 1.0)
これをPoisson方程式の係数行列を作る際の境界ウェイトとして使用する。境界速度はSWEグリッドからバイリニア補間して取得する。
・Neumann条件：
grad_p_dot_n = rho * dot(u_star - U_boundary, n) / dt_apic

Step 6: 圧力Poisson解法とG2P
・Jacobi-PCG法により3次元の圧力を解く。Fractional Volume phi を考慮した係数行列を用いる。
・laplacian(p) = rho/dt_apic * divergence(u_star)
・収束条件：||r|| / ||b|| < 1e-4、最大反復 80、G2P実行。
・グリッドから粒子へ速度と3x3アフィン行列を書き戻す。

6. --- レンダリング統合 ---
・SWE層: 水深 h をハイトマップ(ディスプレイスメントマップ)としてテッセレーションメッシュを生成する。
・APIC層: Screen Space Fluid Rendering (SSFR) で深度から法線を再構築する。
・視覚的統合: 両者の法線をブレンドし、境界には砕波フォーム（泡）を重ねる。

7. --- 数値安定化方針 ---
・乾燥セル（h < h_min）は計算をスキップし、負水深はHLLCフラックス内で厳密にクランプする。
・孤立粒子は非圧縮計算の不安定源となるため、フェードアウトして削除する。

8. --- 設計思想の整理 ---
・厳密保存は質量のみ。
・運動量は可能な限り引き継ぎ、エネルギーは視覚的な「勢い」として扱う。
・全体が一つの流体として動いているように見せることを最優先する。

9. --- 注意事項 ---
・難解な手法から逃げて簡単なアルゴリズムで誤魔化さず、本仕様書通り実装すること。
・値のクランプや不自然な減衰といった場当たり的な対処（ハック）は避け、アルゴリズムレベルでの根本的な解決を原則とする。
・仕様の変更はよく吟味する。特に、保存則に反するヒューリスティックな手法の導入には慎重になること。
・Unityでは内部でcbufferとして扱われるため、コード自体に書く必要はない。
起こりうるバグの多くは圧力が機能しない、粒子が暴れるといったものです。これをそのままAIに投げるとあちこちで値をクランプするハックを提案されますが、これは本質的な解決になりません。実際、原因の多くはUnityの最新のURPの情報をAIが持っていないことや、小難しいミスに囚われてもっと手前の大きなミスを見逃していることでした。

そのため、原因の切り分けは自分で行い、AIにはその後の修正のみを任せるほうが結局時間の節約になると感じました。

 課題と展望iGPUでは描画のほうがボトルネックになってしまい、これ以上突き詰めても得られるものが少ないと判断したため、現状で完成としました。もっとよいパソコンが手に入ったら修正したい、現段階での課題を載せておきます。
SWEの範囲が狭く、FFTとの境目がバレやすい
地形のポリゴンの隙間から水が少し見えていることがある
地形とマウス以外のインタラクションが未実装
水際の泡の表現が上手くいかないことがある
さらに、今回はバージョンアップごとに仕様が大きく異なっていて上手く作ることができませんでしたが、SSFRでの描画に切り替えることで圧倒的な高速化が達成できると思います。
	FFT	SWE	MAC	SPH	PIC	FLIP	APIC
インタラクション	×	〇	〇	〇	〇	〇	〇
波の巻き込み表現	×	×	△	〇	〇	△	〇
安定性	◎	〇	◎	△	〇	△	〇
数値粘性の低さ	粘性なし	△	△	△	×	〇	〇
GPU計算効率	◎	〇	△	×	△	△	FLIPの約1.5倍
特徴	触れない	1次元低い	格子法	粒子法	格子＋粒子	PIC改良	FLIP改良
カテゴリ	処理名	所要時間
Compute	SWEシミュレーション	1.3ms
	APIC 圧力計算	5.2ms
	ボクセル化	29.2ms
	その他(移流,P2G,G2Pなど)	45.1ms
Render	不透明オブジェクト	4.6ms
	透明オブジェクト	40.3ms
	その他(SSAO等)	16.3ms
合計		106.3ms
手法	最小処理時間(※1)	最大処理時間(※2)	備考
フルAPIC	約7,700ms(※3)	約8,900ms(※3)	FFT領域まで全てAPICで計算した場合の推定
APIC+FFT	約168ms	約186ms	昇華・凝縮のシステムがない場合の負荷
REIMU(本作)	約83ms	約106ms	APIC-SWE-FFTカップリング
SWE+FFT	約78ms	約91ms	APICの立体的な波を諦めた場合の負荷
フルFFT	約47ms	約47ms	インタラクションを諦めた場合の負荷
4. 既存の流体シミュレーション手法の比較と限界

各手法の特徴と評価

MAC (Marker and Cell)

SPH (Smoothed Particle Hydrodynamics)

PIC (Particle-in-Cell)

FLIP (Fluid-Implicit-Particle)

APIC (Affine Particle-in-Cell)

SWE (浅水方程式 / Shallow Water Equations)

FFT (高速フーリエ変換 / Fast Fourier Transform)

手法まとめ

5.手法の選定

各手法の構成と採用アルゴリズム

APICの内部構成

SWEの内部構成

FFTの内部構成

6. 異なる流体を繋ぐ遷移の実装(Coupling)

通常のLODからの脱却

SWE → APIC（昇華：Sublimation）

昇華判定

鉛直方向への運動量変換

APIC → SWE（凝縮：Condensation）

凝縮判定

運動量の加算

保存則と安定化の工夫

SWE ↔ FFT間の遷移

7.Compute Shaderのアーキテクチャとデータフロー

ファイル構成とデータ・命令の流れ

ディスパッチ（Dispatch）のフロー

8.性能

検証環境

設定

動作速度

フレームタイム内訳

既存のソルバとの比較

9.技術的課題と解決策(FAQ)

10. おわりに

AIを活用したモダンな開発フロー

課題と展望

Discussion
