UE4のNiagaraのSimulationStageで簡易な物理シミュレーションを行う

#概要
簡易な物理シミュレーションの実装を通して、NiagaraのSimulationStageの使い方を共有します。
題材には、3体問題で有名な3体の質点系の万有引力シミュレーションを選びました。

結果として、以下のような動きが実装できました。
![3BodiesSimulationNiagara.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/0e114e92-6602-784c-6a6f-56035cfc73ca.gif)
地味すぎるエフェクトですが、そもそも絵としてよいものを出すことを目的とはしていないので、そこはご了承ください。:blush:

#背景
NiagaraにSimulationStageが追加されると聞いたときに僕がまっさきに思ったのは、「手軽にコンピュートシェーダで物理シミュレーションができるのか！　すげぇー！」ということでした。
そのためSimulationStageを触りはじめたのですが、自分の結論としては、物理シミュレーションに十分使えるという感触を得ました。

今回SimulationStageを触る過程でいろいろな資料にあたりましたが、情報がまだ少なく、上記のような簡単な結果でさえ出すまでにいろいろと試行錯誤しました。
SimulationStageはUE4.25からExperimentalとして入ってきた若い機能なので、初心者向けの情報はまだ少ないように思います。
また、公開されているサンプルはビジュアル的に凝ったものが多く、すぐ物理シミュレーションがやりたい人にとっては情報過多に感じました。
そのため、今回自分がやったような至極簡単な物理シミュレーションは、コンパクトなチュートリアルとして共有の価値があると思いました。
（まあすぐにでもEpicさんから優れたチュートリアルが出てくるのではないかと思いますが。:blush:）

僕は今回SimulationStageを触り出す前は、Niagaraはチュートリアルやって一部のソースコード読んだ程度の状態でした。
ですので、Niagaraに慣れた人からすると変な実装をしている部分があるかもしれません。
そのような個所を見つけたら是非ご報告ください。

#前提知識
UE4のドキュメントにあるチュートリアルをこなした程度のNiagaraへの理解があるものとします。
以下ページにある、操作ガイドのところにあるチュートリアル達です。
[Niagaraチュートリアル達](https://docs.unrealengine.com/ja/RenderingAndGraphics/Niagara/index.html "Niagaraチュートリアル達")

#SimulationStageとは何か
SimulationStageは、Niagaraに設けられた新しいステージです。
EmitterSettingsステージのEmitterPropertiesモジュールのSim TargetでGPUCompute Simを選択していれば、「Add Simualtion Stage」ボタンでいくつでもSimulationStageを追加できます。
![EmitterProperties.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/fa1b7edb-937f-b1ff-22b4-1ebd1d33c2ef.png)

下の画像のUniversalGravitationという名前のステージは僕が追加したSimulationStageです。
GenericSimulationStageSettingsモジュールでは、IterationSourceでは何に対して処理を行うかを指定し、Iterationsではサブステップ数を指定します。
![NE_MultiBody.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/942dbf3d-b8ec-a313-2b23-a1793aaaa381.png)

SimulationStageにあって既存のParticleUpdateステージにない特徴は、以下だと思っています。

1. Iterationsでサブステップ数を指定できる
1. 複数のSimulationStageを追加できる
1. IterationSourceで処理する対象を選択できる

上から順番に説明します。

既存ステージであるParticle Updateでは1フレームに一回の処理しかできません。
物理シミュレーションでは1フレームに一回のシミュレーション計算だと精度的に不足することが多く、サブステップ分割できる機能は必須です。

複数のSimulationStageを追加できるのは、多段階の処理を重ねていくことが多い物理シミュレーションでは必要な機能です。

IterationSourceではParticlesとData Interfaceというものを選択でき、後者によってParticle Updateステージでは不可能だったシミュレーションアルゴリズムの実装が可能になります。
Particlesを選択した場合は、Particle Updateステージと同様に、処理単位がパーティクルごととなります。
Data Interfaceについては、今回のサンプルでは扱わないので、この記事では解説しません。

#プロジェクトとソースコード
[Githubレポジトリ](https://github.com/monguri/UE4NiagaraSandbox "Githubレポジトリ")
CPUシミュレーションをするためのソースコードはSourceフォルダに入れています。
Content/MultiBodySimulationに今回作ったレベルとNiagaraアセットを入れています。

以降の節ではソースコードとNiagaraアセットを順に解説していきます。

#CPUシミュレーション
いよいよ3体の質点系の万有引力シミュレーションを実装します。
いきなりNiagaraでGPUシミュレーションで実装してもデバッグがしにくいので、正常動作しているかどうか簡単に比較できるようにCPUシミュレーションを先に実装します。

結果として、以下のような挙動を実装しました。
見栄えがするような初期位置、初速、万有引力定数を与えています。
60FPSで、1フレームを10サブステップで分割しています。
![3BodiesSimulation.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/35d8121c-a52e-cf3e-e5e2-6e3f1c7e7ffb.gif)

万有引力は距離の逆二乗で表せます。
[Wikipedia](https://ja.wikipedia.org/wiki/%E4%B8%87%E6%9C%89%E5%BC%95%E5%8A%9B "Wikipedia")

シミュレーション処理は、MultiBody.cppのTick関数内にあります。
抜粋したものを貼っておきます。

```cpp
// サブステップの時間
float Delta = DeltaSeconds / NumIteration;

for (int32 IterCount = 0; IterCount < NumIteration; ++IterCount)
{
	for (int32 i = 0; i < MassPoints.Num(); ++i)
	{
		AMassPoint* MassPoint = Cast<AMassPoint>(MassPoints[i]);

		for (int j = 0; j < MassPoints.Num(); ++j)
		{
			if (i == j)
			{
				continue;
			}

			AMassPoint* AnotherMassPoint = Cast<AMassPoint>(MassPoints[j]);

			// 位置変数の前サブステップの値を用いて速度を計算する
			const FVector& Diff = AnotherMassPoint->Position - MassPoint->Position;
			float DistSquared = Diff.SizeSquared();
			const FVector& Acceleration
				= Gravity * AnotherMassPoint->Mass
				/ FMath::Max(DistSquared, 1.0f) // 0除算にならないよう、最低距離を1cmとする
				* Diff.GetSafeNormal();

			MassPoint->Velocity += Acceleration * Delta;
		}
	}

	// 速度から今サブステップの位置を更新する
	for (int32 i = 0; i < MassPoints.Num(); ++i)
	{
		AMassPoint* MassPoint = Cast<AMassPoint>(MassPoints[i]);
		MassPoint->Position += MassPoint->Velocity * Delta;
	}
}

// アクタの位置を更新する
for (int32 i = 0; i < MassPoints.Num(); ++i)
{
	AMassPoint* MassPoint = Cast<AMassPoint>(MassPoints[i]);
	MassPoint->SetActorLocation(MassPoint->Position);
}
```

#NiagaraのSimulationStageを使った実装

次にNiagaraで同じものを実装します。
CPU実装の処理に該当する部分を先に見せたいので、今回は、内部シミュレーション処理を先に説明し、次にそれを実行するためのNiagaraアセットの設定を説明する、という順番とします。

以下が内部シミュレーション処理のHLSLです。

```cpp
OutPosition = InPosition;
OutVelocity = InVelocity;

if (InDeltaTime < 1.0e-4f)
{
	return;
}

#if GPU_SIMULATION
bool bValid;
float3 Position;
InAttributeReader.GetVectorByIndex<Attribute="TmpPosition">(InParticleIndex, bValid, Position);

for (int i = 0; i < InTotalSpawnedParticles; ++i)
{
	if (i == InParticleIndex)
	{
		continue;
	}

	float3 AnotherPosition;
	InAttributeReader.GetVectorByIndex<Attribute="TmpPosition">(i, bValid, AnotherPosition);
	float AnotherMass;
	InAttributeReader.GetFloatByIndex<Attribute="Mass">(i, bValid, AnotherMass);

	float3 Diff = AnotherPosition - Position;
	float DistSquared = dot(Diff, Diff);
	float3 Acceleration = InGravity * AnotherMass / max(DistSquared, 1.0f) * normalize(Diff);

	OutVelocity += Acceleration * InDeltaTime;
}

OutPosition += OutVelocity * InDeltaTime;
#endif // GPU_SIMULATION
```

これをカスタムHLSLノードで実行しています。
カスタムHLSLノードでないとパーティクル同士の相互作用計算をループでやることはできないのではないかと思っていますが、もし通常ノードだけでやる方法をご存じの方がいればぜひ教えてください！

おおよそ、CPUシミュレーションとロジックが一致しているのが確認できると思います。
細かな書き方については、ContentsExampleのNiagaraAdvancedレベルの3.6 Position Based Dynamicsのサンプルをおおいに参考にさせていただきました。

forループの階層がCPUシミュレーションのときと違って一階層だけになっていますが、これは、このカスタムHLSLノード自体が各パーティクルに対して実行されるのと、サブステップのループはSimulationStageで自動で行われるからです。

InAttributeReaderにはParticleAttributeReaderを渡しています。
ParticleAttributeReaderは、あるパーティクルの処理の中で他のパーティクルの変数値を取得するための機能です。
パーティクル同士の相互作用を実装するときには必須の機能となります。
他のエミッタから出たパーティクルであってもアクセスできます。
SimulationStageだけでなく、ParticleUpdateステージでも使えます。

一回のサブステップで、全パーティクルの位置の更新を同期させたいので、前サブステップの位置をキャッシュしておくためのTmpPositionおよびTmpVelocityという変数を各パーティクルにもたせています。
位置の更新を同期したいのは、同じサブステップ内で、パーティクルによっては先に位置更新されたものが別のパーティクルの計算に使われたということを避けたいからです。
CPU実装では、内側階層のfor文を2つに分割することで対応していました。

（もしかしてNiagaraが内部で別に一時変数をもっていてサブステップごとの更新を同期してくれてたりしないかな？と思ったのですが、今回はそこまで調査できていません。
わかる方がいたらぜひ教えてください！）

次に、SimulationStageでこのHLSLを実行するためのNiagaraアセットの設定を順を追って説明していきます。

Niagara Systemはこんな感じになっています。
![NS_MultiBody.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/9a23e2f8-98e0-38b5-7ecc-ab932db0b912.png)

UserParametersモジュールで、Userネームスペースに追加した変数のデフォルト値を設定しています。
Userネームスペースに追加した変数の値はレベル上に配置したNiagaraActorから設定することもできます。
![UserParameters.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/be308e00-9472-2043-eb4a-a0148a33ef3b.png)

後からカスタムHLSLで使うために、Paticle Attribute ReaderをEmitterUpdateステージで作っておく必要があります。
![SetParticleAttributeReader.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/95aaebbb-cc66-e7a6-706e-4e65b7856940.png)
Paticle Attribute Readerは複数作れるため、Emitter Nameがユニークなキーになっています。

ParticleSpawnステージのInitializeParticleでPosition、MassをUserネームスペース変数から直接指定しています。
![InitializeParticle.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/30f642a4-6d3c-eebb-ae62-d493b8e76515.png)
Niagaraでは初期位置などの初期値を決めるときは範囲ランダムや何かしらのソース位置を使うことが多いですが、今回は再現性のある物理シミュレーションを行いたいので入力値を直接設定する必要があります。
それには、上の図のように、ModeはDirect Setを選択し、Select Xxx from Array、インデックス指定はDirectでReturn Exec Indexを選択します。

配列から設定する場合はReturn Exec Indexを使うことで配列のインデックスとパーティクルのインデックスを対応付けられます。
Particle SpawnステージにおいてExecIndexとはパーティクルのインデックスと思ってもらっていいです。
これは、Particle Updateステージでも同様ですし、Iteration SourceをParticlesにした場合のSimulation Stageでも同様です。

初期速度の設定はAddVelocityで行っています。
色についてはSelect Xxx from Arrayを使える機能が既存モジュールになかったので、Set Particle Colorsというモジュールを自作して設定しています。

Particle Updateには、Solve Forces and Velocitiesモジュールを追加していません。
AddVelocityモジュールを追加すると「Solve Forces and Velocitiesモジュールがいるのではないか？」という感じの文章の警告が出て、Fix Issueボタンを押すとSolve Forces and Velocitiesモジュールは自動で作られます。
今回はDissmiss Issueボタンを押して、あえて警告を無視しています。
これは、今回は位置と速度の更新をSimulation Stageでのみ行いたいためです。

SimulationStageのGenericSimulationStageSettingsは前述のとおりの設定です。
![NE_MultiBody.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/942dbf3d-b8ec-a313-2b23-a1793aaaa381.png)
Iterationsについては、本来は直値を入れるのでなくUserネームスペース変数を設定したいのですが、方法を見つけられませんでした。
方法を知ってる人がいたらぜひ教えてください！

HLSLの説明のときに説明しましたが、今回は、サブステップごとに全パーティクルの位置更新を同期するため、前サブステップのPositionとVelocityをTmpPositionとTmpVelocityに一旦コピーしています。
![SetTmpPositionVelocity.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/43de0696-1ab2-9708-6651-27a1f82e2813.png)

それを利用して今サブステップでPositionとVelocityを更新するのが自作のGravity Simulation Tmpモジュールです。
![GravitySimulationTmp.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/26de3827-1c37-bda2-6cc4-9a3baa3d4ae9.png)
見ての通り、前述したカスタムHLSLノードを実行しているくらいなものです。

最後にこのシミュレーション結果をMeshRendererに用いて、このNiagara Systemをレベルに配置すると冒頭のgif動画の挙動になります。
![3BodiesSimulationNiagara.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/16365/0e114e92-6602-784c-6a6f-56035cfc73ca.gif)

CPUシミュレーションと見た目の挙動が一致しているので、とりあえずは正常動作しているということにしておきましょう。
ちゃんと比較するなら毎フレームの位置と速度を監視したいですね。
Niagaraで、そのようなデバッグ機能をご存じの方がいたらぜひ教えてください！

以上で、今回作ったものの解説は終わりです。

#参考資料

- UE4.26のContent ExampleのNiagara Advancedレベル

現状、SimulationStageについては資料が少ない状況なので、まずはNiagara Advancedレベルのサンプルを見るのがよいと思っています。
コメントが丁寧に入っているので学習しやすいです。

- [Niagara Advancedを覗いてみる](https://qiita.com/yeczrtu/items/4a846b98e89daa4e5090 "Niagara Advancedを覗いてみる")

SimulationStageについての解説はありませんが、Niagara Advancedレベルの各サンプルについて解説されています。

- [EpicのAsher Zhuさんの動画](https://www.youtube.com/watch?v=-cKgZrrBJ2w "EpicのAsher Zhuさんの動画")

自作の流体のサンプルについて解説しつつ、Simulation Stageの学習素材について語っておられますので参考になります。

#あとがき
コンピュートシェーダによる計算をノードベースのGUIツールでやれるなんて、僕はNiagara以外に知りません。
Niagaraマジリスペクトです。
他に同じようなことができるツールを知ってる人はぜひ教えてください！

今回Niagaraで物理シミュレーションをやったのは、それ自体も目的ですが、NiagaraとSimulationStageの内部挙動を調査する足掛かりにしたかったというのもあります。
これから調査していくつもりです。
気が向いたら、それについてもQiita記事にするかもしれません。

目的が目的なので、今回作ったものは絵としてすごく地味です。
今回の実装を発展させれば、万有引力であれ、PBDであれ、SPHであれ、N体のシミュレーションをやるのは素直にできるのではないか（？）と思います。
興味ある方はぜひアーティスティックな絵を作っていただき、物理シミュレーションエフェクト界隈を盛り上げていただければと思います。:blush:

Nが大きくなった時は、相互作用計算を総当たりでやっていると負荷が大きいです。
そのような用途のために、SimulationStageと同タイミングでNeighborGrid3Dという空間分割機能がNiagaraに導入されています。
ContentsExampleのNiagaraAdvancedレベルで使い方が丁寧に説明されているので、参考にしていただけるといいと思います。