PyMCを使いたくて苦闘。
MCMCを理解する過程で必要となったのですが、これが動きません。
PyMC3 が動いたと教えていただいたライブラリVer.の組み合わせを真似し、仮想環境を構築したのですが、ダメでした。PIP を混ぜてPyMC4にもしてみましたが、ダメ。Anaconda を最新版に入れなおしてもダメ、gcc 関連のパスを通してもダメ。お手上げでした。
よく引っかかったのは Theano。これ、Deep Learning にも流用できる数値計算ライブラリとして有名だったのですが、既に開発が終わって時が経っています。
ちなみに、PyMC4 は TF を利用。TFも２年ほど前はよく利用されていましたが、今は PyTorch に首位の座を奪われているように感じます。
バックエンドが変更になるたびにコードが動かなくなるのは避けたいですね。Docker しかないかな？
Accelerating Fortran DO CONCURRENT with GPUs and the NVIDIA HPC SDK
https://developer.nvidia.com/blog/accelerating-fortran-do-concurrent-with-gpus-and-the-nvidia-hpc-sdk/
昨日のブログです。
OpenACC のデータ転送が面倒で困っていたところ、この記事を見かけました。
All data movement between host memory and GPU device memory is performed implicitly and automatically under the control of CUDA Unified Memory.
これでACCと同程度のパフォーマンスが出たらラッキー。頻繁にメモリへのアクセスがあるので過剰な期待はしていないのですが、いくらか早くなれば嬉しい。早く試したい！
Docker はまだ 20.9 (SDK 単体も)。公開まであと1か月くらいでしょうか？
**************************************
20201217追記
先日、Docker イメージが公開されました。
早速試してみたところ、-stdpar=multicore では OpenMP と同じ結果、速度でした。
が、-stdpar=gpu では発散。うまく計算しません。
コンパイル情報を見てみると、gpu では DO CONCURRENT 内の DO LOOP も自動並列化で GPU に載せられていました。これかな？
OpenACC を利用しようと考えていました。
ACC の実装は比較的手軽だそうですが、場合によっては CUDA に並ぶことがあるとのこと（以前、チャレンジしたときは失敗しました）。
私はプロではないので、C & CUDA によるゴリゴリチューニングはできません。で、OpenMP からの ACC で、あわよくば狙いです。
OpenACC を gcc がサポートしているようなので使おうかな？などと考えながら、いつもの PGI Community Edition をチェック。以前のVer.と変わらないなあ、と製品版の価格をチェック。すると、以下の表示。
PGI Professional 製品の販売は終了しました。（詳細）( ﾟДﾟ)？！
PGI Compilers＆Tools は NVIDIA ブランドの NVIDIA HPC SDK ソフトウェアとして、新しく生まれ変わりました。NVIDIA HPC SDK
俄然、やる気になりました。使ってみましょう。
**************************************
20201113追記
使ってみました。
dockerイメージを走らせただけ（nvidia-docker を使わなくなっていました）。その上で何の問題もなくコンパイルできました。
計算中に GPU を使っています。kernels 構文を差し込んだだけですので、遅いのですが。
チューニングには時間を取られそうですが、制御は拍子抜けするくらい簡単でした。ありがたいですね。
SPHスキームの長所と短所
圧縮性（SPH法）
長所：陽解法のため逐次代入計算により時間更新を行う。計算が簡単。
短所：粒子の重なりや抜けが生じやすい。結果として数値的不安定性を招きやすい。
→安定させるためにタイムステップを短くする必要がある。（計算負荷増大）
→人口粘性の導入が必要。ただし、適度に設定しないと訛りすぎる。
弱圧縮性（WSPH法）
長所：陽解法。疑似非圧縮性流体の計算が可能。
短所：計算の安定性のために0.5%程度の圧縮性を許容することが必要。
→体積保存性に欠陥
非圧縮性（ISPH法）
・SPH法の微分演算子とMPS法の計算アルゴリズムの組み合わせ。逆も可。半陰解法。
長所：非物理的な人工粘性を導入する必要がなく安定した計算が実行可能。
短所：圧力場の解に激しい攪乱を伴う
→インターバル平均で攪乱は消失。物理ベースCGなどで利用可能。
→防波堤に作用する衝撃波圧など瞬時の値を扱う問題には使えない。
※陽解法、半陰解法のいずれにおいても、安定性向上に有効なのが粒子再配列（高精度粒子法として分類される）
・XSPH
・Particle Shifting：J＝-v∇Ｃ
・Optimized Particle Shifting
弱圧縮性SPH数値スキームでは、圧力場がスプリアス振動を示します。
状態方程式を用いた場合、密度の変化が非常に小さくても、圧力の変化は指数関数的に生じます。それらは、物理的にはありえないような巨大な圧力を生じさせます。
対策
Density filtering: Shepard Filter (0th order), Moving Least Squares (1st order)
人工拡散項の追加：δ-SPH, Fourtakas et al. (2019)
Riemann solver：Godunov SPH
δ-SPH
WSPH法の数値的不安定（引張不安定）の制御手法として、δ-SPH法が近年よく用いられる。δ-SPH法では密度勾配をTaylor級数の一次精度で算定することにより拡散項の評価精度を向上し、自由表面の跳躍を抑制している。
後藤仁志「粒子法 連続体・混相流・粒状体のための計算科学」より
Molteni, D. and A. Colagrossi. (2009)A simple procedure to improve the pressure evaluation in hydrodynamic context using the SPH.
S.Marrone et al.(2011) δ-SPH model for simulating violent impact flows
Fourtakas et al. (2019)
Local uniform stencil (LUST) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models
Godunov SPH
Riemann問題の解を適用し、必要最低限の粘性を自動的に導入する手法。あらたなパラメータ設定が不要。
Inutsuka et al. (2002) Reformulation of Smoothed Particle Hydrodynamics with Riemann Solver
SPH 空間補間が、平滑化長を h とすると、 𝑂(ℎ2) に等しい精度で収束。
1Dの場合、
\begin{align*}
\left\langle A\left(r\right)\right\rangle=\int_{\mathrm{\Omega}}{A\left(r^\prime\right)W\left(r-r^\prime\right)dr^\prime}\end{align*}
A(r’) を r 周りでテーラー展開すると
\begin{align*}
A\left(r^\prime\right)\cong A\left(r\right)+\frac{\partial A\left(r\right)}{\partial r}\left(r-r^\prime\right)+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r^2}\left(r-r^\prime\right)^2+・・・\end{align*}
SPHでは
\begin{align*}
\left\langle A\left(r^\prime\right)\right\rangle=A\left(r\right)\int_{\mathrm{\Omega}}{W\left(r-r^\prime\right)dr^\prime}\end{align*}
\begin{align*}+\frac{\partial A\left(r\right)}{\partial r}\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}
\begin{align*}+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r}\int_{\mathrm{\Omega}}{\left(r-r^\prime\right)^2W\left(r-r^\prime\right)}dr^\prime\end{align*}
\begin{align*}+・・・\end{align*}
1．正規化：影響域で積分すると１になる。
\begin{align*}
\int_{\mathrm{\Omega}}{W\left(r-r^\prime\right)dr^\prime=1}\end{align*}
2．偶関数：中心に対し点対称
\begin{align*}
W\left(r-r^\prime\right)=W\left(r^\prime-r\right)\end{align*}
\begin{align*}
\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}
\begin{align*}=\int_{\mathrm{\Omega}}{\left(r-r^\prime\right)W\left(r^\prime-r\right)dr^\prime}\end{align*}
\begin{align*}=-\int_{\mathrm{\Omega}}\left(r-r^\prime\right)W\left(r-r^\prime\right)dr^\prime\end{align*}
\begin{align*}=0\end{align*}
3．\begin{align*}q=\frac{r-r^\prime}{h}, w\left(r-r^\prime\right)=\frac{\alpha}{h}f\left(q\right), 𝑓: ℝ+ とすると\end{align*}
\begin{align*}\alpha\int_{\mathrm{\Omega}_q} f\left(q\right)dq=1\end{align*}
\begin{align*}dr=dq・h\end{align*}
\begin{align*}
\alpha\int_{\mathrm{\Omega}_q}{\left(r-r^\prime\right)^2f\left(q\right)dq}=\alpha h^2\int_{\mathrm{\Omega}_q}{q^2f\left(q\right)dq}\end{align*}
1, 2, 3より、
\begin{align*}
\left\langle A\left(r^\prime\right)\right\rangle=A\left(r\right)+\frac{1}{2}\frac{\partial^2A\left(r\right)}{\partial r}\alpha h^2\int_{\mathrm{\Omega}_q}{q^2f\left(q\right)dq}+・・・\end{align*}
\begin{align*}=A\left(r\right)+O\left(h^2\right)\end{align*}
カーネル関数の基本的な性質
