2021年10月31日日曜日

PCを使う?

先日、以下のようなやり取りがありました。

私の質問

  • 新しい共有PCを2か月くらい使用しようと考えている。
  • 1日1回程度、計算状況を確認するだけで、席は空いている。
  • そちらの浸透流解析で何コア使うのか?メモリーはどの程度使う予定か?影響しないように開けておく。

他部署の先輩の回答

  • 浸透流用に導入予定のソフトが非力な PC では動かないケースを懸念している。
  • 新しいPC はハイスペックなのでそれを利用する予定だった。それが使えなくなると非常に困る。

先輩にとっては難しい質問だったようです。知らないことには答えない、という典型的な回答でした。
どうも浸透流は初めてのようで、支配方程式や計算負荷を御存じない様子。そのためか計算はプロに任せ、結果をもらってソフトで表示して、報告書にまとめるだけのようでした。それなら、コアは1個でメモリも数㎇でしょう。ソフトの要件が「メモリ8GB以上」でしたのでどのような PC でも動くでしょうし、「1GB以上のビデオメモリを実装したOpen GL 4.5対応のグラフィックカード」ですので、古い Quadro で問題ないと思われます(新しいPC は GeForce です)。
まったく心配ないのですが、初めての知らないことで御心配だったのでしょう。2,3回やり取りしましたが「使ってくれるな」を繰り返されるだけでしたので、先輩ということで理詰めを控え、引きました(若い方からすると、こういうのを老害と思われるのでしょうか)。

ヒトは誰しも知らないことについて構えてしまうものですが、過大に備えて周りに迷惑をかけることは避けたいものです。今回、それを避けるには PC やソフトを知ることでした。普段から使用している道具の知識を update しておけば、周りに迷惑をかけることはなかったと思われます。無知は罪です。

道具に使われることのないよう、他山の石として肝に銘じておきましょう。

**************************************
ちなみに、このソフトの浸透流ソルバーは私が部分改変してコンパイルしています。Win11になればコンパイルしなおすべきなのですが、御存じないのでしょうね。ま、計算は御自分でされないみたいなので問題があること自体に気付かれなくても支障ないのですが。

20220123追記
ハイスペックPCを確認すると32コア。動かししたいと言われていたソフトは並列化されていません。計算されないものの、ソルバーも並列化非対応。
情報を仕入れず、知識をアップデートしないまま仕事をされると、迷惑なのは周りとお客様です。無知は罪です。

波の伝播

振動の伝播を解こうとして試行錯誤しています。

最初は簡単なモデルで試そうと考え、数十㎞の矩形領域の表面に点震源を与えました。

まずは CTC さんの SoilPlus でモデル化。しかし、計算しだすと完走までおよそ数か月かかることが判明。

これはダメだということで、Misdas さんの FEA NX に切り替え。ソルバーが GPU 対応です(なぜか GeForce 推奨)。
モデルを作って境界条件を張る際に気づいたのですが、側方の粘性境界が実装されていませんでした。伝達境界を張ろうとしましたが、コチラは3次元に対応していないようで晴れません。モデルを作る前に確認しておけばよかった。

で、FLAC 。
こちらは、node に 荷重を与えることが難しいようです。Fish で組み始めましたが、未だ完成せず。

最後にOpenSWPC。
これは以前、MPI が引っかかりコンパイルできなかったのですが、Ver5.20 になっていたためか今回はすんなりと通りました。サンプルを動かして動作確認OK。
実体力(物体力)を時系列で入力するのは左右対称のモーメント時間関数を利用し、点震源を重ね合わせることで対応できそうです。2D の例題で試してみました。

 
Boxcarで10秒間一定の力をかけた場合

 
Cosineで2秒毎(1秒間のラップ)の重ね合わせで10秒間

一定の実体力(Boxcar)をCosineで再現したのですが、内側の違い(9~10秒の関数形の相違)を除けば、一致しているるように見えます。振幅の比較は sac ファイルを確認ですね。ま、最終的にはコードの確認が必要でしょう。


こういった波の伝播に関しては、地震分野のプロがいくつかのコードを公開されていました。が、いずれも震源は断層をモデル化する方式で、時系列の力の入力はサポートされていません(当前なのですが)。OpenSWPC における重ね合わせの利用も、想定されていないでしょうね。

今のところ、FLAC が遅ければ OpenSWPC しか選択肢が残りません。これが正しく機能しないとわかったら、次に行きましょう。

2021年10月24日日曜日

Kratos

セミナーで 紹介されていた Kratos を試してみました。
https://github.com/KratosMultiphysics/Kratos

GiD を1か月フリーで使用できるとのことでしたので、それを利用。
Kratos に限らず、OpenFOAM なども扱えるようです。一通りのジオメトリファイルをインポート可能で、メッシュも nastran などに対応しているようでした。

MPM に関しては、Example に2D が含まれていました。
メッシュを切って run するだけでしたが、?でした。メッシュを切った時点で、解析領域と円の領域の接点が一致していません。円の方が粒子の発生させる領域なのかな?と思っていましたが、メッシュを切っていますし、よくわかりません。

うーん、どういうことなのでしょう。
うーん、寝かせておきましょう。時間が解決してくれそうな気がします。


Unless Making Progress

今日、セミナー動画を見ていると MPM が出てきました。
国内では流行っていないなあ、と思いながら検索。すぐに清水さん、大林さんの文献が引っ掛かりました。掲載が2012年ですから、意外と古い。私自身も2014年に書き残していますね。https://phreeqc.blogspot.com/2014/11/mpm-material-point-method.html


先日、他社さんの資料を見ていると、プレッシャーメーターからcφを求めている方がいらっしゃいました。これ、基準書に乗っているのですが、実務で報告さているのを初めて見ました。
https://phreeqc.blogspot.com/2011/01/blog-post_8311.html

 

前者の流行らなかった一因は PC スペックだろうと思われます。が、それも今は改善しつつあります。
後者は思考の停滞でしょう。組織慣性が働いているのかもしれません。

先達が示された通り、現状維持は後退と同義です。変化しないと進歩しません。
次の10年、新たなことに挑戦し続ける人が増えたら良いと期待します。

2021年10月22日金曜日

AMD + ifort

他支店で、AMD の32コアを積んだ PC を購入したとのこと。

まだ GPU に対応しておらず CPU 並列化のみのソフトは多いので、一つの選択肢かもしれません。が、OS が Windows10 のみ。以前コンパイルしたソフトを動かそうとしたが動かず困っているという連絡でした。

AMD Ryzen は手元に実機がありません。AMD + UBUNTU での並列化はよく聞くのですが、Winはナシ。手探りで始めることに。

エラー内容を見ると、どうやら Intel CPU用の命令が邪魔をしているようでしたので、まずそれを外すことに。そのままでは変更するオプションが多いので、コマンド打ちで0からオプションを追加する方針にしました。

何度か試行して、スタック領域を拡張し忘れていたことに気づき、追加。これで動きました。

今回、ifort ではコンパイルできたのですが、gfortran ではダメでした。何が違うのかわかりませんが、前者の速度と安定感は以前より感じています。

このあたり、いつかプロに御教授願いたいところです。

 

2021年10月20日水曜日

Animated Graph

 
SPH の動画を見ていて、ふと、グラフをどのように動かしているのか?と疑問に。

 


何で作られたのだろうと調べると、Matplotlib にその機能がありました(灯台下暗し)。

from matplotlib import animation
interval=100/3 #300 frames * interval 100/3ms = 10s

ani = animation.ArtistAnimation(fig, frames, interval,repeat=False)
ani.save('./Force.mp4',writer="ffmpeg",dpi=300)

for loop でグラフの作画を進めた後(frames リストに格納した後)、 animation.ArtistAnimation で集約。ani.save でmp4 として保存しました。

SPH の計算結果を動画にするのは blender。と思っていたのですが、これが落ちまくり。png で書き出してから動画にすることにしました。
apng だと PowerPointが対応していないため、これもm4へ。cv2.VideoWriter でOK。30 fps で10秒の動画とし、グラフと合わせます。

できた2つの mp4 を PowerPointに貼り付けたら、同時にスタート。
できました。


 

2021年10月17日日曜日

サンプリング周波数

空間周波数について調べていたのですが、サンプリング周波数として欲しい内容が書かれていました。空間周波数も同様ですのでメモしておきます。

https://www.keyence.co.jp/ss/products/recorder/lab/voltage/point.jsp

サンプリング周波数÷信号波形周波数 記録波形への影響
10倍以上 正確な波形表示・記録が可能
2倍以上10倍未満 波高値が小さくなります。波形形状が荒くなります。
2倍未満 エイリアシングが発生し、実際の波形とまったく異なる波形を表示。

関連→https://phreeqc.blogspot.com/2020/10/2_13.html

GEORAMA でボクセルメッシュ

GEORAMA v3.3
ソルバーの方の GEORAMA です。使用するのは初めてです。

複数の csv 標高データを SoilPlus へ取り込もうとしました。が、そのままでは薄層が多すぎてメッシングが困難になることが目に見えています。
で、事前に Civil3D+GRORAMA でモデリングしようとしたのですが、なぜかLandXML を境界面として選択できませんでした。
結果、csv 読み込みが可能な GEORAMA 単体を選択。

メッシュの大きさに比べ薄い層が多くある場合、ソリッドやメッシュを切ることは難しいのが一般的です。この場合、割り切ってボクセル出力を選択せざるを得ません(解析の種類にもよりますが)。GEORAMA ならボクセルの重心付近の地質をあてはめたメッシュを自動で書き出してくれます。
また、異なる層が同標高となっている場合、他のソフトではある程度ダミーの層厚を持たせる配慮をしないと次の解析に持っていけませんでした。GEORAMA のボクセルなら、これも気にしなくて済みます。割り切ると便利です。

750,000メッシュの作成に1時間くらいでしょうか。その後は nastran 形式で作成し、Soilplus で読み込み(これも時間がかかります)。

無事にメッシュを作成できました。

***********************
20211019追記

ボクセルメッシュは重心位置の地質が適用される仕様だそうです(サポートさん談)。メッシュ内での各地質の体積は考慮されていません。残念です。


2021年10月16日土曜日

DesignSPHysics

セミナーで紹介されていた DesignSPHysics。
https://design.sphysics.org/

以前より公開されていることを知っていましたが、使っていませんでした。
今回、その機能を体験したのですが、これがなかなか素晴らしい。FreeCAD に 組み込むマクロ(Python)ですが、よくできていました。基本的なシミュレーションなら、これを通してジオメトリの取り込み、計算、可視化を容易に実施できます。

復習を兼ねて、手元にある3Dモデルに水を流してみました。ジオメトリは Civil3D で作成していましたので、まずはこれを取り込みます。

変換手順
1.Civil3D:ジオメトリ全体を iges 書き出し(部品毎の stl 書き出しは座標が飛ぶので不可)
2.FreeCAD:iges 読み込み → 部品毎に stl 書き出し(座標、スケールを保つ)
※正面図で書き出すこと。書き出したview の上方向がz軸となる
3.FreeCAD+DesignSPHysics:stl 読み込み(「Import GEO」 で1倍で取り込む)

 

読み込んだジオメトリに対し流体か boundary かを指定した後、諸々の計算条件を設定します。
「Execution Parameters」で人工粘性を組み込んだり、計算時間を決めたり。
「Object order」で粒子の発生順序(領域の重なる中で、どちらを後に発生させて上書きさせるか)を決めたり。計算を CPU で実施するか GPU で実施するかも、ここで指定できます。

 パラメータを決めたら、「Save Case」でモデルを保存。
「Run GenCase」の後、「Run」でDualSPHysics を実行します。

現行  Ver. では mDBC 未対応。書き出した xml を直接修正すればOKでした。DEM, Chrono も今後とのこと。これだけ簡単になると使い続けたくなります。楽しみですね。

計算後は「ComputeForce」で壁等にかかる力を抜き出せます。boundary の MK は+10で指定すれば良いみたいですね。どこかに書かれているのかもしれませんが、わからなかったので VTK 書き出し後に MK を確認しました。
VTK書き出しは「PartVTK」、「IsoSurface」など。ParaView のパスを「Setup Plugin」で設定していたら、自動で立ち上げることも可能です。地味に便利な機能です。
VisualSPHysics を使う場合は、Blender でstl(boundary) + VTK(fluid)取り込み。Civil3Dではジオメトリを m単位 で作成していたため、stl を0.001倍指定で読み込むと VTK と座標が一致します。
残念ながら VisualSPHysics は開発が止まっているようで、 2.90LTS 未対応。2.83LTS に落としたら、問題なく動きました。

で、肝心の計算結果は、イマイチ。うーん、残念でした。

粒子法はメッシュを切らなくて済むので、手軽です。DesignSPHysics でより手軽になったかな。今後の開発が楽しみです。あらためて開発チームに感謝!

衝撃力のサンプリング

玉石がコンクリート壁に衝突する際の衝撃力を解析結果から抽出していました。

最初に衝突する6個の石の速度変化と一緒に図化すると、なぜか2個の衝撃力が取れていません。再計算しながら解析画像を書き出し、礫の衝突するタイミングや配置を追認してみましたが、きちんと壁に衝突しています。
1/1,000秒の書き出しではダメ?ということで1/10,000秒で書き出してみました。

はい、取り出せました。小さいのでΔtも小さいのかな。
1/10,000秒ではピーク値が取れていないように見える箇所があり、最終的には1/100,000秒で書き出して評価。この辺り、数値実験では容易です。

もし、実際に衝撃力を計るとすれば、このような高サンプリング周波数対応の応力測定器を探すだけで一苦労でしょう。ないかもしれませんね。

気をつけましょう。

************************************
20211019追記

https://www.jstage.jst.go.jp/article/sabo1973/32/1/32_1_40/_article/-char/ja/
衝撃力の時間に関し、「砂防ダムに対する土石流衝撃力算定とその問題点」で紹介がありました。内部の歪みは50μsec~1msec以内で消えるとのこと。今回は玉石なので、さらに短いのでしょう。1/1000秒間隔だと、拾えない衝撃力があるのも納得です。

2021年10月9日土曜日

Strontium isotopic composition

Strontium isotopic composition of hot spring and mineral spring waters, Japan
https://doi.org/10.1016/0883-2927(91)90053-R

温泉・鉱泉のストロンチウム同位体比が地質により異なるといった内容。
第四紀および新第三紀の火山地域の温泉水の87Sr/86Sr比は0.703-0.708の範囲を示し、花崗岩、堆積岩、変成岩地域のもの(0.706-0.712)よりも低いそうです。

SPH セミナー

昨年、パルマ大の SPH のセミナーを聞いていたのも、ちょうど今頃でした。

今年も開催されています。2週間の予定で、今日で前半戦が終わりました。
相変わらず英語が分かりません。説明は数式ベースなので何とかなるのですが、質疑応答がさっぱり。動画を見直せば何か引っかかるかもしれませんが、英語をクリアしても理解できない内容のような気がします。
参加されている皆さんの理解の速さに圧倒されています。否応なしに、立ち位置を確認させられています。

今年は dx と h 、それに伴う誤差の理解に焦点があてられているようです。これはありがたいですね。今まで、dx は計算速度から決定し、h=1.3dx を当然のように使っていたのですが、エラーに関して気にしたことがありませんでした。求め方すら知りませんでした。以後、気を付けましょう。

理解促進のため、毎日コーディングの宿題が出ています。提出義務はないのですが、比較的容易な内容なので毎日実施しています。手は動かすもので、理解したつもりだった点、あいまいだった点が見つかります。
support domain = khと、カーネルに使用する h=1.2dx の係数をあわせがちになっていたのですが、今回クリアになりました。初歩のガウスカーネルを選択して書いていたのですが、h が分散のところに入るので、カーネルの大きさに直に効くことにあらためて気づきました。セミナーでは、この大きさに誤差が関連することを説明されていました。
近傍探索アルゴリズムは、 以前、Fortran で組まれたコードを読んでいただけでした。今回、何も見ずに組んでみると迷いますね。あの部分、速く処理するためにどうしてたかな?などと。説明では、side effect としてメモリアクセス、キャッシュでの優位性を話されていました。が、さすがにこれを理解されていた受講生は少なかったような気がします。

翌朝の宿題に関する議論では、言語を問わないとあって FORTRAN, Julia, Pythonで書かれた方々が説明されていました。Julia が 出てきたのは面白い。グラフ化が必要ですからね。

あと4日。有意義な時間が続きますように。

2021年10月4日月曜日

df.plot()

Pandas 偉い。

index を datetime しておけば、

date_time rain(mm/h) Q(mm/h)
2009-10-07 01:00:00 0.019231 0.018576
2009-10-07 02:00:00 0.000000 0.018576
2009-10-07 03:00:00 0.019231 0.018576
2009-10-07 04:00:00 0.000000 0.018576
2009-10-07 05:00:00 0.038462 0.018576
2009-10-07 06:00:00 0.019231 0.018576
2009-10-07 07:00:00 0.019231 0.018576
2009-10-07 08:00:00 0.076923 0.018576
2009-10-07 09:00:00 0.057692 0.018576
2009-10-07 10:00:00 0.096154 0.021672
2009-10-07 11:00:00 0.076923 0.027864
2009-10-07 12:00:00 0.019231 0.034056
2009-10-07 13:00:00 0.038462 0.040248
2009-10-07 14:00:00 0.038462 0.043344
2009-10-07 15:00:00 0.038462 0.043344
2009-10-07 16:00:00 0.134615 0.043344

この呪文で
years=[2010,2011,2012,2013,2014,2015,2016,2017]
for year in years:
	df[df.index.year==year].plot()

こうなる。

 

指定した年毎に一気に書いてくれます。 

plot() の中に matplotlib の制御キーワードを渡せます(matplotlib の wrapper)。

Pandas 偉い。 


2021年10月3日日曜日

タンクモデルのパラメータ最適化

先輩&後輩が、SCE-UA でタンクモデルのフィッティングを行っていました。

数年前にSCE-UAを利用した際、その恣意的な合わせこみ方法に引っかかっており、もっと一般的な最適化手法を適用すればどうなるの?と考えていました(手は動かしませんでしたが)。

良い機会ですので、何かコードが転がっていないか?と探し始めたところ、すぐに引っかかりました。

Python 流出解析用タンクモデルのパラメータ検索に挑戦(日流量解析)
https://qiita.com/damyarou/items/63480fb3e49ff2496a76

「パ-ソナル・コンピュ-タのためのタンク・モデル・プログラムとその使い方(第2報)」に沿って作られたようですね。オリジナルに乱数による範囲設定は入っていたかな?

一般的といえるかどうかはわかりませんが、タンクモデル本家の最適化アルゴリズムなので、ひとまず動くかどうか試してみることにしました。

動かないところを直しながら、少し手を加えて完成。動作確認もOK。結果が表示されました。Python だと計算はやや遅いのですが、図化まで一気にできるのが利点です。図化までのtotal 時間を考えると、小さいプログラムならこちらの方が速いので、手放せません。

で、結果はイマイチ。使われていたデータを眺めた段階で「ダメだろう」と思われたので、致し方ない結果。

ま、動きましたので、第2報の内容を読み返してみましょう。


大量処理

マージすると5億行になるデータを  Python で加工しいました。

for loop の速度が遅いことは有名ですが、まずは素直に書いて試行。
が、処理が終わるのに140日ほどかかる見込み。これはダメ。

joblib で parallel にするとともに、アルゴリズムを改良。これで、4時間程度に納まる見込み。
でしたが、メモリ不足を吐きます。SSD 500GB を仮想メモリに使うも足りず、HDD 2TB をたせば、、、フリーズ。頑張れ、Windows!

結局、小分けして処理。
が、これも夜中にフリーズしていました。Windows 嫌い。

最終的には9時間かかりました。計算は速いのですが、書き込みがネックになりますね。

他にも2週間計算していたPCが途中で再起動していたり、1週間かけて作っていた数億個のデータも途中でストップしていたり。

大きな処理は、PC ではなくクラウドが良いのでしょうか?データ転送がネックになるので、悩むところです。プロはどうされているのでしょう? 


ArcGIS + ゼンリン

 ESRI さんが ArcGIS online でゼンリンさんの住宅地図を扱えるようにしたと案内が来ました。

これは嬉しい、と価格を問い合わせました。
結果、現ゼンリンさんとの契約額と1桁以上の差がありました。これは手が出ません。
ゼンリンさん提供の web 経由サービスでなく、数十倍の金額を払っても ArcGIS 経由で利用したい、という業種やニーズが思い浮かびません。123 で背景に使えると便利、という以上のメリットとなると、どのような利用法なのでしょう?

ArcMap が廃止予定となって実質サブスク化したり、サポートのレスポンスが1~2日程度かかるようになったり。ESRIサンには良い印象を持てなくなりつつあった時期に、このサービス提供。まだ優勢的地位は崩れず勢いがあるのでしょう。