Siv3DでSPH(粒子法)の可視化をしてみた

# はじめに

この記事はSiv3D Advent Calendar 2024の24日目の記事です。

みなさんこんにちは。@emadurandal と申します。TypeScriptとC++で自作3Dライブラリを開発している者です。

自作ライブラリ勢として、Siv3Dの存在はかねてより気になっておりまして、今年ついに手を出してみました。
国産の3Dライブラリは数が少ないので、応援とコード勉強も兼ねて触ってみることにしたのです。

中身を見てみると、C++の最新仕様をうまく取り入れた優れたコードに思わず唸りました。
これはぜひ何かに使ってみたい！

悩んだのはネタ（題材）です。何に使おうか。色々考えた結果、SPH(粒子法)に挑戦し、その可視化でSiv3Dを使ってみようと思いました。

# SPHとは

SPH（Smoothed Particle Hydrodynamics）は、流体や気体などの連続体を粒子としてモデル化し、シミュレーションを行う手法です。この手法は、1977年に天体物理学の分野で提案され、現在ではコンピュータシミュレーションでも広く使われています。

## 1. 基本的な考え方
### 粒子による流体表現
流体を小さな粒子の集合として表現します。これらの粒子は物理量（位置、速度、密度、圧力、温度など）を持ち、流体の性質を表現します。

### 近傍粒子との相互作用
各粒子の物理量は、自分の近くにある粒子（近傍粒子）との相互作用を通じて計算されます。この「近傍」とは、一定の範囲（カーネル関数の影響範囲）内に存在する粒子を指します。

## 2. カーネル関数の役割
カーネル関数は、粒子間の影響を滑らかにするための重み関数です。ある粒子の物理量を計算する際、このカーネル関数を使って近傍粒子からの寄与を重み付きで合計します。

## 3. 流体方程式の離散化
SPHでは、流体力学の基本方程式（ナビエ-ストークス方程式など）を粒子間の相互作用に基づいて離散化します。

連続の式（質量保存）: 粒子の密度を近傍粒子の情報から計算。
運動方程式（運動量保存）: 粒子に働く力（圧力、粘性、外力など）を計算し、速度や位置を更新。
エネルギー保存: 必要に応じてエネルギー方程式も計算。
## 4. SPHの特徴
### 長所:

- メッシュ（格子）を使わないため、自由表面や複雑な境界条件に強い。
- 自然に3次元の流体をシミュレーション可能。
- 局所的に粒子を追加・削除できる柔軟性。

### 短所:

- カーネル関数やスムージング長の選択が結果に大きな影響を与える。
- 粒子数が多い場合、計算量が膨大になる。
- 境界条件の扱いが難しい。
- GPGPU化が比較的難しい。

最後のGPGPU化が難しい（特に近傍探索）というのがネックですね。最近CG界隈ではSPHより格子法と粒子法のハイブリッド型で近傍探索が不要なMPM(Material Point Method)の方がメジャーらしいです。挑戦しようとしたけど不勉強すぎてMPMの方は無理でした。許して……。

# スクリーンショット

動いてますね。はい。並列化しているので動きはゆっくりですがリアルタイムで動きます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/12753/db7af759-5382-14b6-08d3-f57c4c166025.png)

https://x.com/emadurandal_cgj/status/1867796229608157695

# 並列化

Parallels STLで並列化しています。std::threadなどと違い、フラグ一つ立てるだけで並列化できるので、ものすごく楽ですね。

```cpp
	std::for_each(std::execution::par, m_particles.begin(), m_particles.end(), [&](auto& p) {
#if USE_LEAPFROG
		Float3 oldAccel = p.accel;
		p.position += p.velocity * m_dt + 0.5f * oldAccel * m_dt * m_dt;
		Float3 newAccel = computeAcceleration(p);
		p.velocity += 0.5f * (oldAccel + newAccel) * m_dt;
		p.accel = newAccel;
#else
		Float3 accel = computeAcceleration(p);
		p.velocity += accel * m_dt;
		p.position += p.velocity * m_dt;
#endif
	});
```

本当はGPUで並列化したかったんですがスキル不足で無理でした許して……っていうかコンピュートシェーダーで近傍探索とかどうすりゃええねん。

と思っていたら、そのものズバリ、コンピュートシェーダーでSPHを実装された記事が登場しました。これみて勉強しなきゃ！

https://zenn.dev/sparkle/articles/217cc2bb44fd9e

# Siv3Dで今回利用した機能

私のSiv3Dスキルはまだ入門レベルですので、本当に単純な機能しか使っていません。
Box、Sphere、Transformer3Dくらいです。

```cpp
void SPHApp::drawField()
{
	Box(m_fieldMin, m_fieldMax - m_fieldMin).drawFrame();
}

void SPHApp::draw()
{

	drawField();

	const Transformer3D transformer{ Mat4x4::Translate(Vec3(-5, -5, -5)) };

	for (const auto& p : m_particles) {
		Sphere(p.position, m_particleRadius).draw(Linear::Palette::Red);
	}
}
```

# Siv3Dの良いな、と思うところ

## 最小限のコードで目的を果たせる設計

C++でアプリ開発なんていいますと、大抵はWin32APIだとかVulkanだとかの話になってきて、main.cppが数百行～千
行超えたりするのもザラだったりしますが、Siv3Dはとにかく最小限のコードで目的を果たせるようになっています。
以下は今回のデモのmain.cppの内容ですが、C++でこんな短いコードに収まるなんて驚きでした。

```cpp
# include <Siv3D.hpp>
# include "SPHApp.h"

void Main()
{
	Window::Resize(1280, 720);
	const ColorF backgroundColor = ColorF{ 0.4, 0.6, 0.8 }.removeSRGBCurve();
	const MSRenderTexture renderTexture{ Scene::Size(), TextureFormat::R8G8B8A8_Unorm_SRGB, HasDepth::Yes };
	DebugCamera3D camera{ renderTexture.size(), 30_deg, Vec3{ 0, 0, -32 } };

	SPHApp app;

	while (System::Update())
	{
		camera.update(2.0);
		Graphics3D::SetCameraTransform(camera);

		// 3D 描画
		{
			constexpr double BoxLength = 10.0;

			const ScopedRenderTarget3D target{ renderTexture.clear(backgroundColor) };

			for (int i = 0; i < 10; ++i) {
				app.update();
			}
			app.draw();
		}

		// 3D シーンを 2D シーンに描画
		{
			Graphics3D::Flush();
			renderTexture.resolve();
			Shader::LinearToScreen(renderTexture);
		}
	}
}
```

## 最新のC++仕様を随所に活用している

最新のC++仕様を活用した実装のため、コードの可読性が非常に高いです。また、C++11/14/17/20の機能で自分があまり使っていなやり方が書かれていたりするので、同じライブラリ開発者として大変参考になっています。

## 導入が簡単

Visual Studioのプロジェクトテンプレートが用意されているので、簡単に始めることができます。
基本的にWindowsで始めるのが一番良さそうなのでしょうかね。
Mac/Linuxだとどういう開発体験になるかも気になるところです。

# Siv3Dで今後ここに期待！

## プリミティブ形状のインスタンス描画に対応してほしい

Sphereなどのプリミティブ形状ですが、今回粒子として大量に表示しているのでインスタンス描画されているといいなーと思ったのですが、まだそこは対応がないようです。

## APIドキュメントが欲しい

APIを調べたくなったのでAPIドキュメントを探したのですが、まだ……ないっぽい？
ただ、ソースコード規模が巷のゲームエンジンみたいに巨大ではないですし、きれいに書かれているのでソースコードを見れば良い、というのはあります。

# 最後に

Siv3Dは素晴らしいライブラリだと思います。

私も最近C++でライブラリを開発し始めましたが、Siv3Dを先輩格としてぜひ色々とお手本にしたいと思っています。

今回は極めて軽い使い方しかできず、期待された方は肩透かしだったかもしれません。申し訳ありません。
Siv3Dはこれからも触って徐々に慣れていきたいと思いますので、どうぞよろしくお願いいたします。

みなさま、良いクリスマスを！