Turbidity current in Monterey Canyon

昨日の文献で示されていた、avalanche/tutorials/montereycanyon の Allrun を読んでみました。

Rauter, M. et al.(2018) Applications of the Finite Area Method on a Geographic Scale: From Dense Snow Avalanches to Turbidity Currents

The second model is based on the work of Parker et al. (1986) for turbidity currents. The solver models a turbulent mixture of sediment and water, including entrainment of sediment at the bottom and water at the top. It will be applied to a turbidity current in Monterey Canyon, off the coast of California.

参照先のソースを継ぎ接ぎし、要点のみ一読できる形にまとめ、見通し良くしました。
まず、この例では海底地形の STL ではなく、Esri ASCII ラスター形式読み込んでいます。ArcGIS で xyzポイント→ラスター変換→ASC保存 をお考えなのでしょう。文献では "GIS data is translated to OpenFOAM dictionaries using python scripts. "とありましたが、その script は見当たりませんでした。どこかにあるのでしょう。
解析範囲はshpを読み込んでいます。

gridToSTL
    system/gridToSTLDict
        //Output file
        stlname "constant/surface.stl";
        //Topography file
        gridname "constant/gisdata/bathymetry.asc";
            ncols        1194
            nrows        620
            xllcorner    574999.06
            yllcorner    457937.36
            dx           24.71
            dy           31.27
            NODATA_value  0
             -166.5 -164.13 -158.49 -154.69...
        .//Choose between boundary: fromShape or fromPoints:
        //Read Boundary from Shapefile
        boundary fromShape;
        shapeBoundary "constant/gisdata/aoi_UTM_N10";
        //Division of an edge of the boundary polygon
        divisions 200;
        //Height of the Domain
        domainHeight 1000.0;
        //Offset to reduce coordinate size
        offset (-574000.0 -4057000 0);

次に、polyhedral mesh と finite area mesh の作成。前者で meshDict を読んでメッシュを作成していますが、次の faMesh との連携がわかりません。文献ではcfMeshに関し"A mesh is generated that covers the volume above the terrain of interest; but only the bottom boundary mesh is utilized."とあり、faMesh が底のメッシュのみを取り出して2Dに変換する役割を担っているのかなあと想像しています。
release area は inlet で処理しているとのことでしたが、その詳細も見つけられませんでした。

pMesh: cfMesh-polyhedral mesh
    system/meshDict
        maxCellSize 200;
        surfaceFile "constant/surface.stl";
                additionalRefinementLevels 2;
                "wall0"
                    newName "inlet";
                "wall(9|10|11)"
                    newName "outlets";
                "top"
                    newName "top";
                "terrain"
                    newName "terrain";
makeFaMesh: finite area mesh
    constant/faMesh/faMeshDefinition

あとは通常運転。
０フォルダの作成（複製)。
restore0Dir: restore initial fields

ソルバーの設定、計算。
$(getApplication)
    system/controlDict
        faParkerFukushimaFoam
        #include "shapefileWrite":
            offset (-574000.0 -4057000 0);
        #include "gridfileWrite"
            ncols, nrows...offset
        #include "autoAreaToVolumeMapping"
            prefix  "fa_";

いくつか不明点は残りましたが、概ね見えてきました。

Gravitational mass flows + OpenFOAM

Rauter, M. et al.(2018) Applications of the Finite Area Method on a Geographic Scale: From Dense Snow Avalanches to Turbidity Currents
https://www.itn-slate.eu/wp-content/uploads/2019/10/ESR_14_MatthiasRauter-Abstract-7thESIOpenFOAM-Berlin2019.pdf

河川・砂防分野でよく使われる２次元浅水方程式にスカラー輸送、沈殿等を加味し OpenFOAM2006 に導入されています。$FOAM_TUTORIALS/incompressible/shallowWaterFoam/ のような簡単な例をベースに、雪崩や土石流などの Gravitational mass flows (GMF) に特化した形で実装されたのでしょう。

手元のWin11＋ OpenFOAM 2012 だと（当然ですが）コンパイルできませんでした。ダメもとでUbuntu20.04 ＋OpenFOAM９でも実施。こちらも Foundation 版には含まれていないソルバーを見に行ったみたいでダメでした。

正攻法の Ubuntu20.04 ＋ OpenFOAM2112 だとコンパイルが通りました。で、tutorial の Avalanche を Allrun。完走しました。
paraFoam で結果を見るとそこそこ。GIS上に持ち込めるようですが、中を読まないとわからないですね。説明が省略されていますので、地形の導入から結果の export まで、中を見たらわかるということなのでしょう。単純な GMF なら、手元の Fortran プログラムや iRIC の方が良いのかもしれません。

OpenFOAM ユーザーが多いのでしょう。多くのソルバーや解析例が公開されています。ありがたいのですが、初心者には敷居の高い open source です。

ObsPy その3

Obspy 1.2.2 で　import　できなくなりました。numpy で引っかかっています。

仮想環境を再構築してもダメ。
1.2.1 に douwngrade してもダメ。

Git で issue を検索するとヒット。
Import error for obspy v1.2.2 (windows, pip install) #2940
https://github.com/obspy/obspy/issues/2940

$ conda search numpy
$ conda install numpy=1.21.5

numpy を落とすと動きました。
不用意に upgrade してしまい、1時間ほど無駄にしました。

SciKit GStat

解析値と点群から作った地形グリッドの差分を可視化し、チェックする必要が出てきました。

数年前は Surfer で対応していましたし、Civil3Dでもよく実施していた内容です。
ただ、解析値が出るたびに手作業が発生するのは非効率。これ以外の計算や図化はプログラムに載せていましたので、この部分もその流れに含めたいところです。

さすがにクリギングぐらいは Python で容易にできる時代になっているでしょう、と思って検索したら、ありました。

SciKit GStat
https://scikit-gstat.readthedocs.io/en/latest/

地球統計学ライブラリは複数ありましたが、今回はこれを選択。必要な機能を最初に確認できたライブラリで、BESTということではありません。
解析値に対し Gaussian model で variogram を計算し、クリギングによる補間、グリッド値を算出します。これだけです。

あとは他のグリッド状の解析値と合成し、LPグリッドとの差分をとり、図化・出力。これで手作業から脱却できました。便利な時代です。ありがたい。

FloPy その2

自身も FloPy を動かしていなかったことに気づきました。

ひとまず tutorialと example をいくつか動かしてみました。
MODFLOW 関連ソフトは含まれているのかと思いきや、別途DLが必要でした（当然）。MT3DMS はUSGS版でない方が必要でした。

動かしてみて気づいたのは、MT3DMS にいくつかの吸着式や離散化スキームが実装されている点。昔、Visual-MODFLOWを触っていた時にも選択できたのでしょうが、気づきませんでした。
結果を見ると、風上差分が訛っています。これ、OpenFOAM の1次精度風上差分でも同様でした。2次精度では改善していましたので、MT3DMS の upstream は前者なのでしょう。MOC（とHMOC）だと訛らない、MMOCだと訛るというのは意外。やってみないとわかりません。

先の揚水量と透水係数の変化に応じた水面変化は、tutorialケースを修正することで容易に確認できました。無理な条件を設定すると収束しませんし、取水深度を考慮する必要があります。水位ではなく水頭コンターが表示されますので、同じ平面位置でも、どの深度に観測孔を設けるとどのような観測水位になるかも想像しやすい。EXCEL よりはるかに多くの情報を与えてくれます。

やはり、その時代に即したツールを使う方が効率的です。

技術者の階段 その2

数年前にお世話になった御客様から電話を頂きました。

どうも、御一緒した仕事の続きをされるみたいでした。
当時、地質統計学を用いた3次元確率分布の可視化を含め、好評だったと記憶しています。覚えていてくださったのはありがたい話です。で、今の担当部署にバトンタッチ。
結果、executive が「こんなことはできないよ」と一蹴されていました。担当部署では技術的に備えておらず、一般的な調査のみで仕事を片づけられたようです。私も当時使っていたソフトは部署移動時に手放してしまいましたのでサービスできず。残念です。

別部署の後輩さんとの雑談で、水位と透水係数の話が出てきました。

水位を変えるために透水係数をどう設定するかが分からないようでした。先輩さんの専門が異なるので、地下水に関する十分な指導を受けていないのでしょう。
昔 EXCEL で作っていた揚水時の2次元水位変化の book を渡しておきました。単純化しているので正確ではないのですが、透水係数や揚水量に応じてどのような水面が得られるかイメージを掴む程度のことは可能でしょう。
本来なら基礎理論を理解し、水理公式や FloPy などでイメージを掴んでから浸透流計算に移るべきです。が、後輩さんの部署は基礎理論はおろかプログラミングも計算もされないようですので EXCEL ぐらいしか選択肢がありません。残念です。

いずれも単純な技術力不足なのですが、それで通用するのがヒトの世。どこかで技術者の階段を必死に登ることをやめたヒトたちが集まって、バイアスがかかり「常識」が生まれ、技術力不足を感じなくなったのでしょう。怖いですね。

しかも、このような状況に陥らないようにするために技術士の総監があったはずです。専門技術のみならず、管理技術を勉強される方が増えると、お客様にとっても喜ばしいことです。そう思えば日々研鑽される方が出てくるのではないでしょか。

一人でも多く、そのような「技術者」と呼べる方が増えることを願います。

CfdOF その2

CfdOFに関しては、多孔質媒体を含む流体解析が目的。

OpenFOAM では Darcy Forchheimer（ダルシー・フォークハイマー）model として実装されています。
https://openfoamwiki.net/index.php/DarcyForchheimer

日本語の文献もありました。　
高田ほか（2019）防波堤基礎捨石より流入する津波や高潮に関する水理実験および数値解析
吉岡ほか（2010）高透水性多孔質体中の非ダルシー流れに関する考察

係数は Ergun の推定式が良いとされていますが、透水係数からでも求められます。今回は一つ目の文献と同条件で多孔質媒体の透水性を定め、テストしてみました。

VOF（Volume of Fluid）法での2相流です。越流からのエアのトラップも表現されていますが、メッシュをもっと細かくしないとダメですね。

粒径を半分にした場合がコチラ。

透気係数や空気の圧縮性が気になりますが、設定を見つけられませんでした。表面張力の設定は transportProperties の中にありました(sigma)。後からエディタで修正が必要でしょう。

OpenFOAM の基本は tutorial をコピーして利用するスタイルかと思いますが、こちらの方が便利。自由表面の流れと地下水を同時に扱う場合、これをベースに設定を詰めていくスタイルの方が楽そうです。

*************************
20220205追記

式を展開していくと、高田ほか（2019）の(12)式とは微妙に異なる形になりました。数値は一致したので、誤植でしょうね。

CfdOF

Surface water と Groundwater にかかわる文献を読んでいて、CfdOFにたどり着きました。

Integral Flow Modelling Approach for SurfaceWater-Groundwater Interactions along a Rippled Streambed
A finite-volume solver for two-fluid flow in heterogeneous porous media based on OpenFOAM

CfdOF は FreeCAD のワークベンチ。
早速、Win10 ＋ FreeCAD v0.19　に CfdOF をインストールしようとしたのですが、以下のエラーが発生。

エラーが発生したためこのワークベンチをインストールすることができません. 先ずは、見当たらないコンポーネントをインストールしてください。

Gitのレポジトリを見に行くと、インストール手順が書かれていました。
CfdOF: A Computational fluid dynamics (CFD) workbench for FreeCAD
https://gitlab.com/opensimproject/cfdof/-/tree/master

Before installing CfdOF, the Plot workbench must first be
installed into FreeCAD using the Addon manager:

Run FreeCAD
Select Tools | Addon manager ...
Select Plot in the list of workbenches, and click "Install/update"
Restart FreeCAD
Repeat the above for the "CfdOF" workbench
For installation of dependencies, see below

このあと、dependencies をチェックする流れ。その中の一つに OpenFOAM が必要でした。これは FreeCAD に含まれているものと勘違いしていました。

Ubuntu20.04 では 既に OpenFOAM を走らせて結果を可視化するまで環境を整えていましたので、そちらへ行移行。したのですが、FreeCAD v0.18 から addon manager を立ち上げても CfdOF が出てこない。調べてみると、v0.18 の不具合でした。

*********ここから読み飛ばしてOK**************
仕方ないのでWin10 に戻り続行。
管理者権限で FreeCAD を立ち上げてから 編集‐設定 で CfdOF を選択。必要なソフトをボタンポチポチでインストールしました。OpenFOAM は CFD 版の2012がデフォでした。Ubuntu 側は Foundation 版でしたので、ちょうど良いかも。

ところが、dependency cheker がうまく走りません。無視して計算を流しても、エラーになります。
因果関係はわからず、これは重そう、ということで再び Ubuntu へ。

*********ここまで読み飛ばしてOK***************
ターミナルから
$ sudo add-apt-repository ppa:freecad-maintainers/freecad-stable
$ sudo apt update
$ sudo apt upgrade

保留されてしまったので、
$ sudo apt dist-upgrade
これで v0.19 になりました。

次は gmesh のインストール。
$ sudo apt install gmsh
ダメでしたので、公式HPからパッケージを DL して /opt の下位へ展開。

OpenFOAM の directory を指定。
/opt/openfoam9
paraview,gmesh,output も同様に指定

cfMesh, HiSA をポチっとインストール（OpenFOAM の directory を先に指定しないと入りません）。

Win版で作成した fcstd データを Ubntu 側にコピーして、実行。
gmsh のケースがないといわれたので作成して再実行。

走りました。これで環境は整いました。

*********ここから読み飛ばしてOK**************
Win10側でも mesh を作成しなおすと走りました。他のケースでも試してみましたが、mesh作成は計算の直前に実施した方が良さそうです。順番があるのかもしれません。
checker がうまく走らない点の原因はわかりませんが、Win10側でも計算環境が整いました。

plotly in Jupyter

Jupyter Lab にて plotly の図が表示されませんでした。

公式HPを見てみると、以下の案内がありました。
https://plotly.com/python/getting-started/

$ conda install "ipywidgets>=7.6"
$ conda install -c plotly jupyter-dash
$ conda install python-kaleido

はい、これで無事に表示されました。

以前は Surfer で作成していた地形変化の図化（色塗りコンター）。今は Plotly でも対応可。xyz の DataFrame にしてしまえば、容易に図化できます。個人的には メッシュ配列に加工しなくて良いため、matplotlib よりも好みです。

import plotly.graph_objects as go
import plotly.offline as offline

fig = go.Figure(data =go.Contour(z=df.z,x=df.x,y=df.y,))
fig.update_layout(title='contou',
                  autosize=False, width=500, height=500,
                  xaxis=dict(ticks='outside'),
                  yaxis=dict(ticks='outside'),
                  margin=dict(l=20, r=20, b=20, t=20))
fig.update_traces(zmin=-100,zmax=100,colorscale= 'Edge')
#fig.show()
offline.iplot(fig, filename='contour', image='png')
fig.write_image('./contour.png')

pyc to py

pyc ファイルを py へ decompile。

使用したのは uncompyle6 v3.8。Python3.10 だと動かなかったので、3.8 で環境構築。これで動きました。 

数10ファイルのうち、エラーはひとつ。他はほぼ完ぺきに再現できているでしょう。

import uncompyle6
import glob

for f in glob.iglob("./**/*.pyc", recursive=True):
    print(f)
    pyc_path = f.split('.')[1]
    py_path = '.'+pyc_path + '.py'
    with open(py_path, "w") as fileobj:
        try:
            uncompyle6.decompile_file(f, fileobj)
        except:
            print("Error")

********************************
20220108追記
再現内容を見るとインデントのおかしな箇所があり、そのままでは完走しませんでした。修正すればOKでしたが、完ぺきではないようです。ま、ここまで再現できれば素晴らしいのですが。windowsだったからかな？

SeisNoise

Timothy Clements, Marine A. Denolle
SeisNoise.jl: Ambient Seismic Noise Cross Correlation on the CPU and GPU in Julia
Seismological Research Letters (2021) 92 (1): 517–527.
https://doi.org/10.1785/0220200192
https://docs.juliahub.com/SeisNoise
https://github.com/tclements/SeisNoise.jl

ambient seismic noise cross correlation のライブラリです。これまで見てきた passive 手法の中では最も新しく、シンプルに書けそうです。しかも CUDA 対応済み。
https://phreeqc.blogspot.com/2021/05/passive.html
https://phreeqc.blogspot.com/2021/06/blog-post_25.html

ちょうど手元の Win10 ＋ Jupyter がエラーを吐いて動かない状態でしたので、再構築がてらインストールしてみました。
まずは、既存の Julia と miniconda をアンインストール。関連フォルダをすべて削除。レジストリもクリーンにしてから再起動。

Julia 1.65 をインストール後、
pkg> add IJulia
pkg> add SeisNoise
pkg> add SeisIO; build; precompile
pkg> add Plots
 
miniconda 3 をインストール後、
$ conda install jupyter lab
$ jupyter lab

規定の browser だと表示されなかったので、起動 browser を変更。
$ jupyter lab --generate-config
~/.jupyter/jupyter_lab_config.py
 c.ServerApp.browser = '"C:\\Program Files (x86)\\Microsoft\\Edge\\Application\\msedge.exe\" %s'

これで立ち上がりました。
Github の README.md のソースをコピペして run。
同じ絵が出ました。ＯＫです。

F-net の sac が手元になかったため、FDSN(IRIS) にある気象庁の静岡と岐阜のデータに変更。周波数を低めに設定し、run。

１日間と短いのですが、そこそこの形になりました。

DL に時間がかかりますが、仕方がないでしょう。
テストは CPU ＋1日間のデータのみでしたが、感覚的に速い方だと思います。GPU を使えるなら、大きなデータを扱えそうです（データを用意するのに時間がかかりますが）。

exsample は今後整備される予定のようです。楽しみですね。

Mahalanobis distance

マハラノビス距離を利用する異常検知手法は、機械学習系の図書に載っています。

これまで、ユークリッド距離での異常検知しか実施したことがありませんでした。が、それらの結果は散々。今回、淡い期待を込めてマハラノビス距離も潰しておくことに。

ガウス分布の指数部分のみなので、計算自体は簡単です。
pandas の DataFrame の形で多次元データを成形してしまえば、分散共分散行列、逆行列、平均を減ずる操作が容易。機械学習系のフレームワークを使用せずとも、pandas と numpy のみで完結します。

で、結果は散々。ま、不可を確認できましたので良しとしておきましょう。