弾塑性SPH その2

新たに頂いたコードを参考に、以下の組み合わせで計算ができるように改造しました。

・弾性
・弾塑性
・弾性‐流体
・弾塑性‐流体
・流体

これらのほとんどは日本語の図書でも紹介されています。個体と流体を同じように実装できるのはSPHならではでしょう。

弾塑性（一部流体）の挙動はコチラ。端から土圧により崩壊し始めます。

少しc・φをあげると、崩れにくくなります。その場合、まず角が楔のように崩れ、次に長辺の2/3が崩れ、最後に残りが崩れるような時間差が見られました。ちなみに、前回のアップしているイメージは「弾性‐流体」の計算結果で、全体がペシャっと潰れるような挙動を示しています。これだけ自由な表現ができると、多様な自然現象を再現することが可能でしょう（予測は相変わらず困難ですが）。

まだ、いくつか疑問点が残っています。資料を探して早く完成させましょう。

壁粒子の圧力

SPHで壁際の粒子が壁に張り付く現象を、後回しにしていました。
ある程度動くようになったので取り組んでみることに。

まずは mDBC で様子見。
効果はありますが、残念ながら完全ではありません。いくつかの粒子は壁に残り、計算時間も余計にかかりました。

負圧が効いているのは想像がつくのですが、他の修正方法を思い付きません。
そこで googling。
すると、「まさにコレ！」といった回答を見つけました。条件も全く同じです。
Trying my chances here; anyone knows why particles stick to the wall in my SPH simulation? - Specific Domains / Modelling & Simulations - Julia Programming Language (julialang.org)

Here the simplest fix is simply to put the boundary particle density to reference density, if it ever gets below that reference value.
In this way a wall can never produce suction, but only “push”.

当たり前すぎて一般的には解説されていないのでしょうね。初心者には思いつかないノウハウです。
密度でなく圧力を変更。数行の追加で計算時間もかからず、簡単に解決できました。

弾塑性 SPH

2020年の秋から始めた弾塑性SPHコードの修正が完了しました。

基本コードを頂いてから寝かせていたため実質1年程度ですが、SPHを知ってからは14年目です。長い。
何とかSPHの離散化や各種補正方法、境界条件について基本的なしくみを理解したと言える証になるでしょう。いや、よかった。 

いくつかのハマり処があったのですが、式の実装時の静水圧Pと静水圧応力σの符号の整理が効いたと思われます。多くの文献を引用する際に、これらの書き方の流儀を統一してから実装に落とし込む必要があったのでしょう。

粒子が凝集する現象を hourglass,  zero energy mode だと勘違いしたり、他にも似たような勘違いで多くの寄り道をしましたが、振り返るとそれらも良い勉強になったと思います。
https://phreeqc.blogspot.com/2022/05/blog-post.html

早くVerification を済ませたいですね。ベンチマークは何が良いのでしょうか？
それが済んだら高速化。倍精度を多く使っていますので、GPUではなくスパコン利用（MPI並列）のほうが現実的でしょうね。

まだまだ先は長い。

相互相関関数の計算

西田究のホームページ (u-tokyo.ac.jp)
解説 (u-tokyo.ac.jp)

地震動に関する様々なシミュレーション結果をお手軽に見ることができる、素晴らしいサイトです。特に、地震波干渉法で相互相関関数が出来上がっていく過程を見たことがなかったので、非常に興味深く見ていました。

　　初期

数十分後

相互相関関数を求める際は、時間領域（直接法）か周波数領域（FFT）を利用します。配列の要素数が多ければ周波数領域での計算が早いでしょう。Python の scipy.signal.correlate では手法を指定できますし、指定しなくとも早い方を自動で選択してくれます（Ver.1.13.0）。 numpy.correlate は直接法のみです（Ver.1.26）。

以下、計算例です。設定次第なのですが、以下だと答えは①≠②=③
①直接法

sig1 = np.std(data1)
sig2 = np.std(data2)
data1 = data1-data1.mean()
data2 = data2-data2.mean()
cross_correlation = np.correlate(data1, data2, mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(data1)
f2 = np.fft.fft(data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum)
                            ) /sig1/sig2/N

③FFT2

f1 = np.fft.rfft(data1)
f2 = np.fft.rfft(data2)
cross_spectrum2 = f1 * np.conj(f2)
cross_correlation = np.fft.irfft(cross_spectrum2
                                ) /sig1/sig2/N

0をパディング + シフトで①＝②となります。

①直接法

padded_data1 = np.pad(data1, (0, N), 'constant')
padded_data2 = np.pad(data2, (0, N), 'constant')
cross_correlation = np.correlate(padded_data1,
　　　　　　　　　　　　　　　　 padded_data2,
　　　　　　　　　　　　　　　　 mode='same'
                                ) /sig1/sig2/N

②FFT1

f1 = np.fft.fft(padded_data1)
f2 = np.fft.fft(padded_data2)
cross_spectrum = f1 * np.conj(f2)
cross_correlation = np.real(np.fft.ifft(cross_spectrum))
cross_correlation = np.fft.fftshift(cross_correlation
                                   ) /sig1/sig2/N