コルモゴロフ・スミルノフ検定

2群のデータを比較するためにヒストグラムを作成しました。
それを眺めながら、分布形状が似ている、似ていないを判断したいのですが、主観だと説明性に欠けます。ということで統計的検定の出番です。

以前、2標本（対応ナシ）の条件で Mann-Whitney's U test を使用しました。
https://phreeqc.blogspot.com/2020/12/mann-whitneys-u-test.html

今回は2群のサンプル数が異なるとともに、サンプルが非常に多く何らかの確率分布を仮定できそうでした。
調べてみると、ありますね。
・2標本のコルモゴロフ・スミルノフ検定
two-sample Kolmogorov-Smirnov(K-S) test

手順は、①2群を同じ階級幅で区分して累積度数を求める。②累積度数を累積確率に変換し、その差の絶対値を計算。③絶対値の最大値を使って検定統計量を求め、④棄却限界値と比較し判定する。棄却限界値は自由度ｆ=2のχ2分布から求めます。

これは、好都合。ヒストグラムを全面積で正規化し、縦軸を確率密度で記載していたので処理は半分終わっています。あとは累積確率に直して差の最大値を求めるだけのようなもの。χ2分布の適用も、おかしくはない分布ですから、これで決定です。というか、2群のヒストグラムを書く際に、ルーチンワークとして求めておけば良い内容です。Python コードを残しておきましょう。

def hist_plot(col,x1,x2,label,bins):
    fig = plt.figure(figsize=(15,5))
    ax = fig.add_subplot(1,1,1)

    freq,bin_list,patch=ax.hist([x1,x2], density=True, alpha=0.5 ,edgecolor='black', bins=bins, label=label)
    ax.set_xlabel(col,fontsize=14)
    ax.set_ylabel('確率密度',fontsize=14)
    ax.grid(which='major',color='gray',linestyle='-')
    ax.grid(which='minor',color='gray',linestyle='--')
    ax.legend(borderpad = 0.8,framealpha=1.0,shadow=True)
    
    #2標本コルモゴロフ・スミルノフ検定（χ２分布、各標本数40以上、対応の無い2標本）
    df_p=pd.DataFrame({col:bin_list[1:],label[0]+'_確率密度':freq[0],label[1]+'_確率密度':freq[1]})
    df_p_sum=df_p.cumsum()*bin_list[1]
    df_p_sum=df_p_sum.drop(col,axis=1)
    df_p_sum.columns  = [label[0]+'_累積確率', label[1]+'_累積確率']
    df_p_sum['D']=abs(df_p_sum[label[0]+'_累積確率']-df_p_sum[label[1]+'_累積確率'])
    df_p=pd.concat([df_p, df_p_sum], axis=1)

    n1=len(x1)
    n2=len(x2)
    Dmax=df_p['D'].max()
    Dmax2=Dmax*Dmax
    T=4*Dmax2*(n1*n2)/(n1+n2)
    if T<5.99146:
        comment='有意水準5％で、差があるとは言えない。'
    else:
        comment='有意水準5％で、差があると言える。'        
    display(df_freq,'Dmax',Dmax,'T',T,comment)

Dが小さいと、差があるとは言えない（H0を棄却しない）

H0：帰無仮説（差がない）
H1：対立仮説（差がある）

Civil3D2022 サーフェス作成時のポイント数

他支店の方から、Civil3Dでの等高線の作り方を教えてほしいと連絡がありました。

いつもはV-nas Clairで作成されているようですが、非常に遅いので困っているとのこと。
すぐにできると思いましたが、沼にハマりました。

取り込むポイントは4800万点。これから等高線を作ってV-nasで表示したいとのこと。聞けば、その方も社内で偉い人？から投げられた作業のようでした。技術者なら必要な範囲を判断し、指示するはずですが。センスのない話です。

最初にハマったのは、Civil3D でサーフェスを保存できない現象。原因はメモリかと思いましたが余裕があるので異なりました。少しずつ点数を減らして限界を把握しました。結果、今回は4200万点まで保存できました。
調べてみると同様の報告がありました。サポートさんも御存じないようでしたが、以前から存在する現象のようです。現行の2023でも解決していないのかな。

次にハマったのが V-nas Clair の制約。Civil3D で作成した等高線を LandKit でインポートすると、「要素の制限値を超えた【ポリライン】があります。分解もしくは削除されました」というダイアログが出ます。サポートさんに伺うと、1本のポリラインに2^15以上の要素が含まれるとポリラインは分解、スプラインは削除する仕様とのことでした。分解すると起終点を閉じたポリラインは開くので、その間の要素が消えたように見えます。が、Civil の等高線は開いたポリラインのため、影響なし。よかった。その前に TrendPoint のスプライン等高線を読んでいたのですが、一部削除された原因がコレでした。

単純作業だったのですが、ハマりポイントが複数あったため時間がかかりました。ま、各ソフトウェアの限界を知ることができたので、良しとしましょう。

 

Software	ポイントファイル数	ポイント総数	サーフェス作成	サーフェス表示	等高線作成	save-reopen	V-nasClairでの等高線Opren	結果
Civil3D 2022.1.3	RCS:1	48,000,000	×：フィルタなしで15分程度。 作成完了した時点でフリーズ。 https://knowledge.autodesk.com	境界のみ	-	※ reopen時に再度サーフェスを作成しているような挙動。 数分かかる。	kome	×サーフェスも等高線も作成できない。
〃	TXT:1	48,000,000	〇：Civil15分、メモリ20GB	コンター	〇サーフェスから取り出し	×：サーフェスを保存できない。 等高線のみなら可。 https://forums.autodesk.com	〇	〇等高線のみ保存可
〃	TXT:1	48,000,000	〇	境界のみ	〇	×	〇
〃	TXT:2	48,000,000	〇	境界のみ	〇	×	〇
〃	TXT:1	42,000,000	〇	境界のみ	〇	〇	〇	◎
〃	TXT:1	36,000,000	〇	境界のみ	〇	〇	〇
〃	TXT:1	24,000,000	〇	境界のみ	〇	〇	〇
V-nasClair 2021.4	TXT:1	48,000,000	〇点群読み込み2分、TIN作成6分：メモリ18GB ※動かす度に再描画するため遅い。	TIN	×：等高線間隔設定まで6分。 その後一晩以上でも作成されず中止。 メモり100GB	〇：サーフェスreopen8分	-	×遅い。メモリーを使う。 等高線を作成できない
TrendPoint Ver.9（9006）	TXT:1	48,000,000	※時間を要するため中止	-	〇：作成に一晩。メモリ33GB	〇点群、等高線のみ。Save20分、reopen10分以上かかる。	×：スプラインに完全対応していない。AutoCAD等で少しづつポリラインに変換する必要あり。	△遅い。

QGIS + PostgreSQL

ネット注文がメインになったことやコロナ禍ということで、大型書店に通う頻度が低くなりました。久しぶりに通うと、気になる図書がちらほら。この週末は2冊の 図書を読んで過ごしました。

そのうちの1冊がコチラ↓
愛知大学三遠南信地域連携研究センター「地域研究のための空間データ分析入門 -QGISとPostGISを用いて‐」

GIS は「習うより慣れろ」でしたので、この種の図書を読んだことがありません（ArcHydro 関連の図書は購入しましたが）。この図書には PostgreSQL との連携が書かれています。どのような場合に RDB を利用すべきか常々考えていましたので、思わず購入しました。

QGIS から PostgreSQL への接続は、何度か実施しています。が、実務レベルで採用したことはありません。SQL ではなく、Pythonで整形してしまってから GIS へ取り込むか、そのまま Python で表示してしまう流れの方が速かったからでしょう。紹介されていた作業内容も、データ加工の部分では C 等で作成されたライブラリで処理する Python の方が速いのかもしれません。

図書では PostgreSQL に shp を取り込み、PostGIS で空間分析まで実施されていました。このようなことができるとは知りませんでした。初めて PostGIS の役割を理解できたような気がします。Python では空間分析を GIS に任せていました。SQLで書いておけば、データ更新時の繰り返し処理の手間を省けるところが利点でしょうか。

単年で完結する作業ではデータベースを更新する必要がないため、RDB まで必要としません。一方、データを随時更新・追加したり、多くの方が データ入力に携わったりする状況では、RDB の方が良いのでしょう。
図書を読むことで、ようやく答えの１つを得ました。

Power Automate

メールファイルの移動、HPのボタン押しなどを Power Automate で実施。

使ってみて気づいたのは、以下の点。

• リモートだと失敗することが多い（遅延の影響？）
• バックグラウンド処理ができない。（常に画面表示。マウスをさわって別の作業をしているとダメ。当然ですが、実行中は'待機'になります。）
• 無料版だと自動実行（定時実行）ができない。（Loop処理で似たようなことは可能。）

後輩君に聞いてみたら、同じところ（自動実行）で困ったとのこと。どうして解決したのか聞いたら、VBAと組み合わせたらしい。微妙です。

地下水～地表流まで

・層流 : i=au
・層流から乱流への遷移状態、乱流 : i=au+bu²
・乱流 : i=bu²

1. 層流（粒径小、速度低）

Darcy：\begin{align*}-\mathrm{\nabla h}=\frac{1}{k}\boldsymbol{u}=a\boldsymbol{u}\end{align*}　h=p/(ρg)より\begin{align*}-\mathrm{\nabla p}=\rho g\frac{1}{k}\boldsymbol{u}=\rho ga\boldsymbol{u}=\alpha\boldsymbol{u}\end{align*}
Darcy-Kozeny Carman:
\begin{align*}a=\frac{1}{k}=\frac{180\nu\left(1-n_e\right)^2}{g{\mathrm{\Phi}^2n}_e^3d^2}\end{align*}　ν=μ/ρより\begin{align*}\alpha=\frac{\rho g}{k}=\frac{180\mu\left(1-n_e\right)^2}{{\mathrm{\Phi}^2n}_e^3d^2}\end{align*}

2. 層流～乱流への遷移状態、乱流（粒径大、高速）

Forchheimer：
\begin{align*}-\mathrm{\nabla h}=a\boldsymbol{u}+b\left|\boldsymbol{u}\right|\boldsymbol{u}\end{align*}
Ergun(1952)
\begin{align*}a=\frac{1}{k}=\frac{150v\left(1-n_e\right)^2}{gn_e^3d^2}
　b=\frac{1.75\left(1-n_e\right)}{gn_e^3d}\end{align*}

Kadlec and Knight(1996):角ばった粒子に適用。
\begin{align*}a=\frac{1}{k}=\frac{255v\left(1-n_e\right)}{gn_e^{3.7}d^2}
　b=\frac{2\left(1-n_e\right)}{gn_e^3d}\end{align*}

3.乱流

Manning：
\begin{align*}u=\frac{1}{n}R^\frac{2}{3}I^\frac{1}{2}\end{align*}\begin{align*}I=\frac{n^2}{R^\frac{4}{3}}u^2=bu^2\end{align*}
h：水頭（m）
p：ポテンシャル（Pa）
ρ：密度（kg/m3）
k：透水係数（m/s)
u：速度（m/s）
α：係数（s/m）
β：係数（s2/m2）
Φ：球形度(0～1、球=1)
d：有効粒径（m）
μ：粘性係数（Pa・s）
ν：動粘性係数（m2/s）=μ/ρ
ne：間隙率
g：重力加速度(m/s2)
n：マニングの粗度係数(m^−1/3・s)
R：径深(m)
I：勾配

OpenFOAM
https://openfoamwiki.net/index.php/DarcyForchheimer
https://phreeqc.blogspot.com/2022/02/darcy-forchheimer-model.html?m=0
\begin{align*}-\mathrm{\nabla h}=a\boldsymbol{u}+b\left|\boldsymbol{u}\right|\boldsymbol{u}\end{align*}\begin{align*}-\mathrm{\nabla p}=\rho ga\boldsymbol{u}+\rho gb\left|\boldsymbol{u}\right|\boldsymbol{u}=\alpha\boldsymbol{u}+\beta\left|\boldsymbol{u}\right|\boldsymbol{u}=\mu D\boldsymbol{u}+\frac{1}{2}\rho F\left|\boldsymbol{u}\right|\boldsymbol{u}
\end{align*}
Darcy:
\begin{align*}\alpha=\frac{\rho g}{k}=\mu D\end{align*}\begin{align*}D=\frac{\rho g}{k\mu}\end{align*}\begin{align*}F＝0\end{align*}
Darcy-Kozeny Carman:
\begin{align*}\alpha=\frac{\rho g}{k}=\frac{180\mu\left(1-n_e\right)^2}{{\mathrm{\Phi}^2n}_e^3d^2}=\mu D\end{align*}\begin{align*}D=\frac{180\left(1-n_e\right)^2}{{\mathrm{\Phi}^2n}_e^3d^2}\end{align*}\begin{align*}F＝0\end{align*}
Forchheimer-Ergun:
\begin{align*}\alpha=\rho ga=\frac{\rho g}{k}=\frac{150\mu\left(1-n_e\right)^2}{n_e^3d^2}=\mu D\end{align*}\begin{align*}D=\frac{150\left(1-n_e\right)^2}{n_e^3d^2}\end{align*}\begin{align*}\beta=\rho gb=\frac{1.75\rho\left(1-n_e\right)}{n_e^3d}=\frac{1}{2}\rho F\end{align*}\begin{align*}F=\frac{3.5\left(1-n_e\right)}{n_e^3d}\end{align*}

PersianSPHでは、式と係数を選択。ソース内で係数変更可。
OpenFOAMでは、D,F 入力で表現。

PersianSPH その3

最後に、土‐水連成を作ってみましょう。
と古いサンプルに手を付けましたが、使えないコマンドがありました。

・dom.SeepageType
・dom.Time
・dom.KernelType
・dom.VisEq

Version 違いでしょうね。
Source を見てみると、SeepageType はParticle に対して指定するようです。無駄なので、後で修正されたのでしょう。Time は Private から出してやるとコンパイルが通りました。良いのかな？他に影響が出たらその時考えましょう。
後半２つは新しい？サンプルの書き方に修正することで通りました。

 土‐水連成のキーとなる変数は2つのようです。

・SWI
・SeepageType 

SWI:

0 => The seepage force + The bouyant unit weight of soil
1 => The seepage force + The surface erosion(Lift+Drag) + The bouyant unit weight of soil
2 => The seepage force + The pore water pressure from water particles
3 => Zero interaction force

1の文献が見当たりません。ソースは以下の通り。河床変動のようにせん断力を考慮しているのだろうと想像するのですが、何から引っ張ってきた式かは追えませんでした。

double Cd = 24.0*(P2->MuRef/P2->RefDensity)/(P1->d*norm(v)+0.01*h*h) + 2.0;
SFt = (3.0/(4.0*P1->d)*P2->RefDensity*(1.0-P1->n0)*Cd*norm(v)*v) *K;
SFt(1) += (P2->RefDensity*(1.0-P1->n0)*norm(v)*fabs(P2->S-P1->S)) *K;

SeepageType :

0 => Darcy's Law
1 => Darcy's Law & Kozeny–Carman Eq
2 => The Forchheimer Eq & Ergun Coeffs
3 => The Forchheimer Eq & Den Adel Coeffs

Kozeny–Carman だけかと思っていましたが、他にもありますね。https://phreeqc.blogspot.com/2021/01/kozeny-carman-equation.html
deとd15の片方しか指定する箇所がなさそうなので、こちらの文献に従って組まれたのだと思われます。
https://www.sciencedirect.com/science/article/abs/pii/S0266352X17302318
係数α, βはソース内指定です。透水係数に応じてソースを変更する感じでしょうか。

土‐水連成の場合、強度、透水性、SWIなど調整すべき項目が弾塑性に比べて増えます。計算時間は単純に倍。２次元でもパラスタが容易とは言えません。３次元ならなおさら。
ということで試算は２次元を選択。堤防の越流破堤を想定しモデルを作ってみましたが、SWI＝2の The surface erosion(Lift+Drag) が効きすぎ。ここを河床変動式に変更すれば、もっと現実っぽくなるのでしょう。
最終的には流入水の勢いで極端にえぐれるか、越流しても壊れないかの設定になり打ち切ってしまいましたが、現実的な実験値と比較すればそれなりの形状は得られるかもしれません。探してみましょうか。

以上で PersianSPH の試算は終了。弾塑性と浸透＋越流の確認は十分とは言えませんので、今後、もう少し試してみましょう。

PersianSPH その2

次は弾塑性。 

見たところ、スタンダードな SPH による解き方を実装しているようです。弾塑性のフレームワークも標準的。

弾塑性のサンプルがなかったっため、DamBreak の水を土に変更してみました。

2次元では問題なく動きます。

3次元ではバクハツ。
これは SPH のデフォなのでしょうか。ただし、zero-energy mode は出ていません。境界部分の過剰な正圧に起因するだけのようです。優秀ですね。
particle penetration を防ぐための Penalty Parameter が効きすぎていると踏んで、1→0.5に下げてみました。
バクハツは収まりましたが、これはダメ。底面で粒子の沈み込みが発生しています。
底面：0.95、側面：0.7程度が良さそうです。

粘着力を変えると、塊のまま移動する土塊とすべり層が表現されて良い感じ。
拡大すると、

せん断変形が孔内傾斜計の累積変位図のようにきれいに表現されています。引っ張りも安定して解けています。

問題は計算時間。
1秒の計算に8コアで4時間。時間がかかります。OpenMPなので、スパコン利用も不可。実務に持っていくには gpu に載せるか、MPI並列にするか。

続く。