PersianSPH

SPH でも弾塑性を扱いたい、ということで PersianSPH の登場です。

まずは、インストール。HPに従って進めます。
https://www.m2clab.com/persiansph/installation
$ wget "https://github.com/mghkorzani/persiansph/blob/master/persiansph_install.sh"
$ bash persiansph_install.sh

ここで止まります。2年前に更新が止まっているようですので、仕方なし。
GitHub のリポジトリにあるコードを zip で落として展開。直下に persiansph_install.sh があるので、これを使いましょう。

$ bash persiansph_install.sh

Installation is completed.
Close all terminal windows and open a new one to take effect the defined environment variables
Please refer to the tutorial to run a simulation.

persiansph フォルダが作成されます。その中の、CmakeLists.txt を書き換えて目的の cpp ファイルを指定（サンプルから一つを選択）。
その後、新しいフォルダを作り、そこからcmake。

$ mkdir (New Path)
$ cd (New Path)
$ cmake $SPH
$ make

実行ファイルが（New Path）に作成されましたので、これを動かすと計算開始。フォルダ内に output ファイルが書き込まれます。ファイル形式は hdf5 と xmf。表示は Visit or ParaView。hdf5 は読めませんでしたが、xmf は両者とも OK でした（ParaView では XDMF Reader or Xdmf3Reader）。2D だとVisitのほうが見やすいかな。

続きは後日。

**********************
20220603修正
ParaView では XDMF Reader が正解でした。Xdmf3フィルターをかけた場合に正常に表示されない場合が出てきました。
hdf5 と xmf はセットでした。xmfを参照し、hdf5からデータを取り出しているようです。

 

Dtransu V2

 ダイヤさんの Dtransu の Version が 2.0になっていました。

大きな変更点は、並列化対応でしょうか。
OpenMP ですので「ようやく」感は否めません。手元にあった ver.1 改良ソースを見ると、OpenMP で並列化したのが11年前。その他いくつか改良を独自に加えていましたが、それらの8割ほどは今回の Ver.UP で追いつかれました。

Post 処理でも ParaView に対応です。これは個人的に嬉しい。
あとは pre だけですね。2D版 は Gmsh 用の変換プログラムが作られていましたが、3D版は「商用ソフトで」となっていました。残念。

１０年以上、改良成果を無償で公表される取り組みは称賛に値します。次の技術者も育っているのでしょうね。

文献の複写

今日、国会図書館に文献の複写を依頼しました。

まだ、紙での提供のようでした。ペーパレスの時代にメールで送ってもらえないのは、資源も時間も無駄です。

著作権が引っかかっているようなのですが、問題になる箇所を詳しく理解していません。
https://crd.ndl.go.jp/reference/modules/d3ndlcrdentry/index.php?page=ref_view&id=1000177352

海外の雑誌は PDF で購入できます。日本でも PDF に対応したサービスはあるようですが、コピーに配慮が必要なのは紙の時代と同じです。
https://jdream3.com/question/copy-etc.html 

国会図書館をはじめ国内で PDF 配布されるようになるのは、まだまだかかるのでしょうか。

RPA

先輩が新入社員に指示を出していました。

「Python でスクレイピングして、トリガーに引っかかったらメールを出すスクリプトを作って。」

新人さんはあっさり仕上げていました。Python ではなく、Power Automateで。話を聞くと、EXCEL と Outlook を使いたいのでこちらの方が簡単かと思ったそうです。

Power Automate を触ったことはありますが、慣れません。RPA に使われるノーコード・ローコード、どうもしっくりきません。個人的には「トリガーに引っかかったらメールを出す」部分を Python で書いたことがあったため、指示を受けたら何も考えずに Python で書いていたと思います。そこに「もっと簡単な方法はないか？」と考えた新人さんの方が偉い。ノーコード・ローコードの類も自由に扱えないと発想が狭くなります。いえ、頭が固いだけですね。

今年の新人さんたちは Python や機械学習に抵抗がないようです。中には SPH プログラムを Fortran で書いていた方もいました（遅れていなくてヨカッタ）。その上、柔軟な発想もできるようです。

害にならないうちに一線を譲りたいのですが、続く人がいません。彼らが一人前になるまでは、持たないでしょうね。

地質の精度

新年度になってから、BIM/CIM 活用ガイドライン（案）の詳細を見ていませんでした。

4月に共通編令和4年4月版の地質に関する箇所をざっと確認し、 昨年度と大きく変更されていないことがわかったので、後回しにしていました。

改めて確認を進めていたところ、地質の照査の箇所で気になる点がありました。「補間アルゴリズム記録シート確認」です。これは昨年度版にも記載があるのですが、今年度は記載例が載っていました。
第一印象は”面倒”です。これ、10層作れば10層分のシートを要求しているのと同意です。流行らないでしょうね。

シートには最適化誤差も記入されていますが、塩野先生の文献を読んだ方でないと記号の意味すら理解できないでしょう。ちなみに、GEORAMA は最適化誤差を吐き出しません。非対応です。

最適化誤差も大事ではあるのですが、印象としては片手落ち。使用した関数に対しどの程度誤差を許容したかよりも、想定した関数が実際の分布とどの程度異なっているか（どの程度適切に表現しているか）の方が大事です。数値解析で言えば verification だけ示して validation に触れていないガイドラインになっているような印象を持ちました。

といっても、実際の分布との差を示すのは難しいでしょうね。いえ、技術的には地球統計学を使って確率で表現すれば良いのですが、不確かさのみを表現するためのモデリングが別途追加になります。地質分布の不確かさを定量的に取り込める設計法が生まれつつはあるものの、現状ではオーバースペックでしょう。そもそも、一般的には統計に耐えうるほどの数のボーリングを掘るわけでもないですから、相応の重要構造物でしか成り立たないでしょう。（余談ですが、地質リスク学会？があるのに地球統計学を利用している話が聞こえてこないのも、このあたりが関係しているのかもしれません。）

とは言え、使用アルゴリズムとパラメータの 伝達は参考になると思います。前段で作成された地質分布に対し、ボーリングの追加箇所以外の断面をなるべく変えたくない場合などに役立つかも。

私ならシートを簡略化して提出するでしょう。
何気なく掲載されたような1枚の図でしたが、引っかかる図でもありました。

Power BI

久しぶりに EXCEL を触りっぱなしでした。

フォーマットの定まっていない数多くの EXCEL ブックを定型に加工しているだけなのですが、なかなか終わりません。ある程度形が定まっていたら Python で、となるのですが、今回はファイルによってヘッダーまで異なため、手作業です。前処理と思って作業をしています。

その中で、座標値をチェックする必要が出てきました。
最初は ArcGIS にテーブルを読み込ませる形で処理していたのですが、座標を加工するたびにテーブルからポイントを作成するのが面倒になりました。

そこで思いついたのが Power BI。
これが正解でした。EXCEL ブックをソースとして指定しておけば、ボタン一つで地図上の位置を更新してくれます。座標値を更新するたびに、Power BI 側の絵を更新。チェックが容易になりました。測地系を指定できない点に目を瞑れば、簡易な GIS として使えそうです。

また、標準で ArcGIS for Power BI が実装されていました。ArcGIS online の地図を利用可能になるようです。多様なニーズがあるみたいですね。

今日で２/3程度が終了。もう少しです。頑張りましょう。

*************************
20220528 追記

全て終了。5000点を超えたあたりから位置表示がおかしくなりました。表示されるポイント、表示されないポイントが出てきます。再起動すると直る場合があります。量なのか、設定なのか。様子見ですね。

九州地整 3DVR コンテンツ その2

先日、九州地整さん主催の UE4 を用いたハンズオンがありました。

先月試した「川づくり」では UE の操作で止まっていたので、良い機会と思い参加してみました。https://phreeqc.blogspot.com/2022/04/3dvr.html

水を張って、土を掘って、草を植える 3工程です。少ない工程ですので時間内に覚えることができました。
残念ながら iRIC の出番はありませんでした。次に iRIC が出てくるなら、また参加したいと思います。

UE4 での作業は、どちらかというと河川技術者よりも環境技術者向きだと思われます。環境を考えた川づくり案に、河川技術者がアドバイスと計算を追加する形で落ち着く可能性も考えられます。

自由な発想で作成できるのが楽しくなるまで、時間がかからないのではないでしょうか。

shp ファイルの測地系

古い日本測地系の5㎞メッシュが必要になりました。

世界測地系（地理座標系のJGD2000）は手元にあります。shp ファイルです。
ArcGIS で日本測地系（地理座標系のTokyo）に投影変換してしまうと、同じ位置にて座標値が変更されます。座標値を変更せずに場所を動かしたいのですが、方法が分かりません。ESRIさんに質問しようかな、QGISで探そうかな、などと考えながら思い出したのが prj ファイルの書き換え。

JGD2000
GEOGCS["GCS_JGD_2000",DATUM["D_JGD_2000",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Tokyo
GEOGCS["GCS_Tokyo",DATUM["D_Tokyo",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

これで動きました。

SPlisHSPlasH その2

SplishSplash のモデルの作り方がわからないと書き残しています。https://phreeqc.blogspot.com/2020/10/splishsplash.html

これ、多様な解法が実装されています。弾性体の Zero-Energy Mode Control まで実装されており、ゴムのような物体でも安定して計算できています。しかもリアルタイムレンダリング。その結果はゲームの一場面を見せられているようです。

クールタイムを挟んでみると、モデルの作り方はすぐに理解できました（しっかり説明されていました）。json or Pythonで下記のブロックを組み合わせれば良いだけでした。https://splishsplash.readthedocs.io/en/latest/file_format.html

• Configuration
• FluidBlocks
• FluidModels
• Emitters
• RigidBodies
• Fluid parameter block
• Animation fields

では、手詰まっているモデルを持ってきましょう、ということで作業開始。
DualSPHysics 用の stl があったので、それを FreeCAD で obj に変換。それを RigidBodies に wall として、FluidModels にparticle 生成範囲として与えます。簡単ですね。

用意ができたので、いざ計算！

はい、爆発します。
デフォの DFSPH が原因だと気づくまで数時間かかりましたが、圧縮を許容する手法だと安定して解けるようになりました。obj の境界上に粒子を発生させているのが一因のようにも見えましたが、他にも原因がありそうでした。

設定を試行錯誤する中でカーネルは選択できたのですが、smoothing length を設定する箇所が見当たりませんでした。代替えなのか、利用する粒子の最大数を設定する箇所がありました。この説明を見つけられなかったのですが、粒子数を少なくし精度を落とす代わりに計算を速くする工夫をしているのかもしれません（別の用途かもしれませんが）。

計算途中から粒子数が減っているようにも見えますし、流れも悪い。気になる箇所がまだまだ多く残っていますが、実装されている手法が多様で入替が簡単です。いずれ改善策はみつかるでしょう。
（一部の手法を除き）計算が速く扱いので、試行も簡単。仕事につかえるか、というと？ですが、どの手法が良いかの試算には使えます。

もう少し、触ってみましょう。

粒子法と連続体力学の基礎

ここ3週間ほど、SPH の図書や文献を読み返し、コードを触っていました。

昨年度、複数の対策を試みたものの改善できなかったコード。しばらく眠らせていたのですが、以下の図書が2月に発売され、仕事の落ち着いた先月に読み始めたことをきっかけに、調査を再開しました。

浅井光輝「明解 粒子法」

クールダウンの効果なのか、あっさり原因が見えてきました。文献を探して解決案も見つけました。

この解決策で利用されていたのが、基準配置と現在配置。連続体力学の基礎でしたね。これは先の図書でも記載されていました。

物理量	物質表示 （ラグランジュ表示）	空間表示 （オイラー表示）	空間表示 ←物質表示	備考
ベクトル	dX	dx	dx=FdX	F:変形勾配テンソル
応力 （テンソル）	S （第2PK）	σ （コーシー）	σ=J-1FSFT	J:ヤコビアン（体積変化）=det(F)
ひずみ （テンソル）	E （グリーン）	e （アルマンジ）	e=F-TEF-1

この図書、粒子法だけでなく、他の手法も比較として概要が記載されています。思いがけず、良い復習になりました。

解決案は見つけたものの、コーディングが面倒。さらに、私の腕では GPU に乗せても実用的な時間に解けるようにまで持っていけません。また寝かせようかと後ろ向きになっているところです。

*****************************
20220515
文献の補正方法を実装するのは意外と容易で、1日で終わりました。が、結果は少し良くなっただけでした。残念。

20240518
コードのバグは全く別のところにありました。符号を整理して式の実装を見直すことで修正が完了しました。
https://phreeqc.blogspot.com/2024/05/sph.html