Dtransu のGPU化 その3

テスト環境：CPUが高性能、GPUは中（FP64は低）性能でアンバランス。
GPU: NVIDIA RTX 4000 Ada(20GB、CC 8.9)
CPU: AMD EPYC 9754(128 コア)

1. 密度流 (v2 sample3.dtr, KAN3=2 密度連成) 

流れ：GPU + 粒子追跡：OpenMP  のハイブリッドとしました。
流れ系 (setelm ~48% + solpcg ~28%の self-time) が構造的に並列で GPU 理想形。ここを GPU カーネル化して大幅短縮です。

粒子追跡 MOVE1 連鎖 (~20%) はバケット探索の IF/GOTO  分岐が多く、GPU だと並列効率が落ちました。これは OpenMP スレッドに流す形が BESTでした。

2. 表面流＋移流分散 

 一番速かったのは ompacc で、密度流と同じ戦略（流れ：GPU + 粒子：OMP）。全構成中で最大のを記録しました。
surface_flow.f90 (D8 地表水ルーティング) の self-time が ~2.3% しかないため、にここは GPU 化不要、host 常駐で問題なし。残りは密度流と同じハイブリッドという構成でした。

当初、なぜ FP64 が弱い GPU でも速くなったのか不思議でした。AIに聞くと、「典型的な arithmetic intensity（バイトあたりFLOP数）が低い 処理だったから」＋「16tだから」とのこと。この辺りをAIに頼るようではまだまだかな。

• 1要素読むごとの計算量はわずか。演算器はメモリ待ちで遊んでいるので足を引っ張らない。
• GPUは数千〜数万スレッドを同時に走らせ、メモリレイテンシをスレッド切替で隠蔽。結果、自分の~360 GB/sをほぼ飽和できた。
• CPU 側の並列効率が悪く、460 GB/sを使い切れず~16t相当で頭打ち。

純 GPU (acc)  が負けるのは、ハードの性能差が大きなこともありますが、分岐の多い粒子追跡が GPU の thread divergence に弱いのも一因。ここをテコ入れするとGPUの方が速くなるかもしれません。なお、途中でSPH の Box Search 手法を取り入れてみましたが、劇的な効果は見られませんでした。
このような計算をしていると、FP64に強いGPUだとどうなるのか、試してみたくなります。

改変コードの公開は配布元が許可されていません。時代の流れでコーディング問題はほぼ解決したので、いずれ配布元でもGPU対応版やMPI版、表流水連携版を公開されると思います。比較は出るまで待ちましょう。

Dtransu のGPU化 その2

Dtransu Ver.2 が手元にあったので、こちらもGPU化に着手。

テストは付属の密度流です。
もともと、OMP対応済みだったので、GPU化のみでした。が、何度かテストするうちに、再現性に劣ることがわかりました。

128tまで実施しましたが、16tで頭打ち。

v1 でも密度流をかけると、同じように再現性を確保できませんでした。テストするモノです。
密度流で PCG の流れが変わるのと、粒子追跡の OMP 化部分が引っかかっているようでした。最終的には v1 で固定しましたが、OMP＋GPU の速度は v2 よりやや低下しました。
最初は v2 GPU＋OMP でパラスタし、最後に v1 で提出、でしょうか。

Dtransu のGPU化

Dtransu のGPU化を実施。

AIさんがプロファイルをとってくれるので、非常に効率的でした。

RTX4000を使っていたのでFP64性能はイマイチ。それでも、2.5倍速くなりました。次は計算用途も考慮して選択しましょう。

ずっと優先度低で眠らせていた課題。今年の頭に外した途端、解決しました。残りは不連続体＋連続体ですが、SPHのGPU化が先ですね。

優先度中：機械学習のスキル増強
優先度低：流体＋個体（不連続体＋連続体）＋振動
優先度低：Dtransu の MPI/GPU 対応
優先度低：地表流＋地下水＋移流拡散

文献：High-resolution raindrop counting via instantaneous frequency sensing

High-resolution raindrop counting via instantaneous frequency sensing on hydrophobic elastic membranes | PLOS One

AI要約

1. 背景
雨粒の数やサイズを測定する装置であるディスドロメータは、気象観測において重要ですが、従来の装置には課題がありました。高精度な標準機器（JWDなど）は高価で、交流電源や専門知識を必要とするため、広範囲な設置が困難です。
一方で、安価な圧電素子を用いたセンサーは「振幅閾値方式」を採用しており、一度衝撃を検知すると一定時間計算を停止（ロックアウト）するため、時間解像度が低く、連続して衝突する雨粒を見逃したり、小さな雨粒を過小評価したりするという欠点がありました。
本研究は、市販のマイクと疎水性弾性膜を組み合わせ、低コストながらこれらの限界を克服する新しい検知手法を提案しています。

2. 手法
本手法は、振幅ではなく「瞬時周波数」の変化に着目している点が革新的です。

• ハードウェア構成: 一般的なPVCパイプの片端に、安価なプラスチック膜（PETやHDPEなど）を張り、ドラムのような構造にしています。この膜に雨粒が衝突した際の音をマイクで集音します。
• 物理的原理: 雨粒が疎水性弾性膜に衝突すると、膜の固有振動（低周波）とは別に、12–18 kHzの非常に短い高周波過渡振動が発生することを発見しました。この高周波信号は雨粒の衝突に特有のものであり、低周波の背景ノイズと分離しやすい特性があります。
• Acoustic Feature Model (AFM): 音声データから、信号の局所的なピーク（山）と谷の位置を抽出し、ノイズを除去した上で「局所的な周波数推定値」へと変換します。
• 機械学習 (1D-CNN): AFMで得られた特徴量を、1次元畳み込みニューラルネットワーク（CNN）に入力します。このネットワークは、風などのノイズと雨粒の衝突を判別し、衝突の確率を算出します。
• 軽量設計: リソースの限られたマイコンやスマートフォンで動作（TinyML）させることを目的としており、モデルのパラメータ数は1000個未満と極めて軽量です。

3. 結果

• 高い時間解像度: 従来の手法が30ms程度の解像度だったのに対し、本手法は約10ミリ秒の時間解像度を達成しました。
• 精度: 評価実験において、予測されたイベントの80%～90%が実際の衝突時刻から10ms以内に収まっており、高い精度で雨粒を数えられることが示されました。
• ノイズ耐性: 振幅を計算に使用しないため、低周波の振幅ノイズ（風など）に対して頑健であり、高い感度を維持したまま背景ノイズの影響を抑えることに成功しました。
• モデルの有効性: DensenetやResidualブロックを用いたCNNモデルは、基本的なモデルよりもノイズ抑制能力が高いことが確認されました。

4. 考察
本研究の最大の特徴は、材料の物理的な減衰に頼るのではなく、信号処理と機械学習によって時間解像度を向上させた点にあります。

• 利点と応用: 非常に低コストな材料で構築できるため、教育現場や発展途上国での気象観測ネットワーク、短期間のキャンペーン調査などに適しています。
• 現在の課題: 強風による高周波ノイズが誤検知の原因となることがあり、さらなるノイズ除去技術や風防の検討が必要です。また、一度の衝突を二重に数えてしまうアルゴリズム上の課題も指摘されています。
• 今後の展望: 今回は雨粒の「カウント」に焦点を当てていますが、将来的には個々の雨粒の運動エネルギーの推定や、平均的な粒径分布に基づいた降雨強度の算出へと拡張することが期待されています。
• 未解明の物理: なぜ軟らかな疎水性膜への衝突で高周波振動が発生するのかという微視的な物理メカニズムは完全には解明されておらず、さらなる研究の余地を残しています。

塩ビパイプに膜を張って雨音を計測するだけという、これ以上ない実用的なアイデアです。

パイプのサイズや膜の張力は低周波の共振には影響しますが、検知に用いる高周波過渡振動（12–18kHz）には影響しないとのこと。また、 PETやHDPEといった異なる材質でも同様の信号が観察されています。身近な材料かつDIYによる物理的な製作誤差に対しても極めて高い堅牢性を持っている点が素晴らしい。

機械学習モデルも最初からエッジAI向きに作られています。今後の雨量推定がどうなるか、期待を込めて待ちましょう。

コード整備 その3

おまけです。

 ・振動計
200Hzサンプリング以上にするとデータ保存時にバグっていたのですが、AIに見せると問題点を指摘してくれました。浅いコピーをやめて深いコピー＆数を増やすことで問題を回避。500Hzサンプリングでも問題なく保存できるようになりました。

・3次元安定解析
これもAIがバグを発見してくれました。学習量の多いであろうPythonで書いていたためか、割と早く修正が完了しました。正しい部分も怪しそうなところは指摘することもありましたが、そこに注意しておけば、コードの照査にかなり有効です。

・RegionGrow3D
MATLABコードをPythonに変換。近い将来できるようになると思っていましたが、もう既に可能だったとは。
1回目の変換で結果はピクセルベースで86％合致。もう少し詰める必要はありそうですが、なかなか賢い。

全体を通してみると以下のような感触。
◎既存コードを他の言語に変換したり修正したりする。
◎コードの照査
◎WEBサイトの解析
〇新たな実装
〇機械学習の自動化

まだまだ性能が上がっていくのでしょう。楽しみです。

AIがAIを鍛える

GWの間、コア写真から柱状図を作成する機械学習モデルの訓練を、ほぼAIに任せきりにしていました。先日記載したとおり、ほとんど人手を加えず、AIの提案を承認・却下するだけ。最後は、気が向いたときにRQD等のスコアの伸びを確認する程度でした。手をかけなくてもスコアが徐々に伸びていくのは、精神的に楽。連休も終わりましたので、ひとまずここで一区切りとします。

昨年まではモデルをどう構築するか、比較対象としてどれだけのモデルを用意するか、最適化をどう進めるかなどを考えながら自ら組んでいました。これらの作業も、今年からはほぼAIに任せられそうです。相手は人ではありませんから、提案方針に何度ダメ出しをしても気を使う必要はありません。納得がいくまで試行錯誤させることができます。文献収集力や実装の速さは人を上回り、初歩的なミスに目を瞑れば、24時間無休で働いてくれる部下ができたようなものです。

一方で、誰でも機械学習モデルを組める状況になったことで、AI活用そのものを差別化要因とするのは難しくなるでしょう。土木分野でも、画像や動画を使った異常検出や点検、数値データを使った予測といった難しいとまでは言えない作業の自動化があります。これらに投資してきた企業が案件を獲得してきた側面もあります。しかしAIの進化により、誰もが短期間で同等のモデルを構築できる時代になってしまったため、こうした先行投資・利益回収型のビジネスモデルは終焉を迎えつつあるように感じます。短期的には実績がある分、有利な状況が続くとは思いますが。

これからは、より一層アイデアと新規性が求められる時代になるでしょう。発想の豊かな技術者こそが重宝されるはずです。そのためには、基礎力を疎かにしないことが欠かせません。若い方々には、便利なAIを活用しつつ、時間をかけてこそ身につく基礎力を養ってほしいと思います。

文献：A multi-layer SPH method for generic water–soil dynamic coupling problems

A multi-layer SPH method for generic water–soil dynamic coupling problems. Part I: Revisit, theory, and validation - ScienceDirect

AI要約

1. 背景
水と土砂の動的な相互作用は、地滑りや土石流といった自然現象から、防波堤やダムの浸透流といった工学的問題まで、幅広く存在します。これらの現象を数値シミュレーションで扱うには、土砂の大きな変形、水面の自由表面流、そして両相間の複雑な相互作用を同時に扱う必要があります。従来の格子法（有限要素法など）では、大きな変形によるメッシュの歪みや自由表面の追跡が困難でした。一方、SPH法（粒子法）はメッシュフリーであるため、これらの問題に適しています。しかし、従来のマルチレイヤーSPHモデルの多くは、空隙率（porosity）の時間・空間的な変化を正確に考慮していなかったり、流体の体積保存に欠陥があったりするなどの課題がありました。
本研究は、混合体理論に基づき、これらの課題を克服した包括的な3次元数値フレームワークを提案することを目的としています。

2. 手法
本研究で提案された手法の核心は、混合体理論（Mixture Theory）マルチレイヤー（多層）粒子で表現するSPHの実装にあります。
数学的定式化（ADFとIDF）: 以下の2つのモデルが提案されました。

• 仮密度ベース定式化 (ADF): 混合体単位体積あたりの質量（仮密度）を用いる手法。
• 真密度ベース定式化 (IDF): 各相の実際の材料密度（真密度）を用いる手法。

特にIDFは、空隙率の空間的・時間的な変化を厳密に考慮できるよう、本研究で新たに質量保存則が導出されました。

支配方程式と構成則:

• 水相: 弱圧縮性ニュートン流体としてモデル化され、状態方程式（EOS）を用いて圧力を算出します。
• 土砂相: 弾塑性構成則（Drucker-Prager降伏基準）を用い、土砂骨格の変形と有効応力を扱います。
• 相互作用力: 浮力（buoyancy force）と粘性抗力（viscous drag force）を考慮します。抗力には、層流から乱流まで対応可能なErgunの式に基づいた2次形式が採用されています。

数値的な工夫:

• 粒子体積の補正: 流体粒子が土砂領域に出入りする際、空隙率の変化に応じて粒子の体積を適切に調整することで、流体の体積保存を確実にします。
• 補正シェパードフィルタ (CSDF): 粒子の一貫性を確保し、自由表面付近のノイズを低減しつつ滑らかな圧力分布を得るための新しいフィルタ手法が導入されました。
• 境界条件: 壁面における粒子の貫通防止と、自由滑り/ノンスリップ条件を適切に扱うための一般化境界粒子法が実装されました。

3. 結果
提案されたモデルの妥当性は、複数の検証ケースを通じて確認されました。

• 体積保存の検証: 水が土砂領域に流れ込む単純な崩壊試験において、粒子体積を補正する手法（VA）が正確な水深を再現できることを示しました。一方で、補正を行わない従来の手法では、流体体積が保存されず不正確な結果となりました。
• U字管内の浸透流: 空間的に空隙率が変化する土砂を通る水の流れをシミュレーションし、解析解と非常によく一致することを確認しました。また、CSDFの導入により、長時間のシミュレーションでも自由表面の乱れが抑えられることが実証されました。
• 重力下での沈下: 水中に沈められた土砂に重力が加わる動的な過程を扱い、水圧と土砂の有効応力が最終的に理論的な静水圧分布に収束することを確認しました。
• 水中地滑りの再現: 実験スケールの水中地滑りをシミュレートし、土砂の堆積形状や、地滑りによって誘発される波（impulsive wave）の形状が実験データと良好に一致しました。

4. 考察
本研究の結果から、以下の点が考察されています。

• ADFとIDFの比較: 両モデルは数学的には等価ですが、SPHの実装においてはIDFの方が数値的に安定し、精度が高いことが判明しました。これは、ADFが空隙率の変化に対して敏感すぎ、圧力計算にノイズが乗りやすいためです。
• 統一的なフレームワーク: 提案手法は、純粋な水領域、乾燥した土砂領域、およびそれらが混合した飽和領域を、単一の方程式系で統一的に扱える点が強力です。
• 今後の課題: 3次元シミュレーションは計算負荷が高いため、GPUによる加速が不可欠です。また、水中地滑りの実験で見られた土砂の急激な体積膨張（細粒分の拡散など）を完全に再現するには、単なる構成則を超えた、水と土砂のより高度な界面混合モデルが必要であることも示唆されました。

密度を変更しないと、計算が複雑にならず理解が楽。IDFがADFに比べて安定であることも実感しました。ただ、地中で体積を膨張させるというのは初期配置が面倒。私は粒子分割・統合で対応しました。発表が2022年ですから、改良案は出ているのかもしれません。