今までは、
import sympy as sp sp.init_printing()
で、使えたのに何かのアップデートの関係で
import sympy as sp sp.init_printing(use_latex='mathjax')
じゃないと表示されなくなってしまいました。 原因は調査中です。
今までは、
import sympy as sp sp.init_printing()
で、使えたのに何かのアップデートの関係で
import sympy as sp sp.init_printing(use_latex='mathjax')
じゃないと表示されなくなってしまいました。 原因は調査中です。
computational-sediment-hyd.hatenablog.jp
computational-sediment-hyd.hatenablog.jp
ほとんど修正していません。 numbaモジュールをimportすることと
import numpy as np import numba
各defの前に以下ように一行付け加えるだけです。
@numba.jit(nopython=True, parallel=False) def calupre(p,u,v,w,Vdir,Vbound,alphaU,dx,dy,dz,nu,axis):
@numba.jit(nopython=True, parallel=False) def caldp(p,u,v,w,acu,acv,acw,dx,dy,dz):
@numba.jit(nopython=True, parallel=False) def SIMPLE(p,u,v,w,Uslip,Vslip,alphaU,alphaP,dx,dy,dz,nu):
コード全体はこちら
上記の関連記事のとおり、numba.jitが使えない条件は多々ありますが、慣れてくるとほとんど引っかからないです。
私は日頃からnumba.jitを愛用しているため、無意識にnumba.jitが適用できるようにコーディングしている気がします。
前記事から久しぶりに新しいモデルを書いたので楽しかったです。ぜひ回してみて下さい。
最近科学技術計算を書いてなかったのでウォーミングアップにこんな感じの3次元キャビティー流れをSIMPLE法で書いてみました。
大きいものはこちら
単純な3次元のキャビティー流れです。
式の導出は以下を参照
computational-sediment-hyd.hatenablog.jp
def calupre(p,u,v,w,Vdir,Vbound,alphaU,dx,dy,dz,nu,axis): # x: axis = 0, y: axis = 1, z: axis = 2 zero = float(0.0) V0 = np.copy(Vdir) nx, ny, nz = p.shape ac = np.zeros_like(Vdir) axp = np.zeros_like(Vdir) axm = np.zeros_like(Vdir) ayp = np.zeros_like(Vdir) aym = np.zeros_like(Vdir) azp = np.zeros_like(Vdir) azm = np.zeros_like(Vdir) b = np.zeros_like(Vdir) act = np.zeros_like(Vdir) if axis == 0: xs, ys, zs = 1, 0, 0 xmb, xpb, ymb, ypb, zmb, zpb = -1, -1, 0, ny-1, 0, nz-1 elif axis == 1: xs, ys, zs = 0, 1, 0 xmb, xpb, ymb, ypb, zmb, zpb = 0, nx-1, -1, -1, 0, nz-1 elif axis == 2: xs, ys, zs = 0, 0, 1 xmb, xpb, ymb, ypb, zmb, zpb = 0, nx-1, 0, ny-1, -1, -1 for i in range(xs, nx): for j in range(ys, ny): for k in range(zs, nz): c = (i,j,k) if axis == 0: uhxp = 0.5*u[i,j,k] + 0.5*u[i+1,j,k] uhxm = 0.5*u[i,j,k] + 0.5*u[i-1,j,k] vhyp = 0.5*v[i,j+1,k] + 0.5*v[i-1,j+1,k] vhym = 0.5*v[i,j ,k] + 0.5*v[i-1,j ,k] whzp = 0.5*w[i,j,k+1] + 0.5*w[i-1,j,k+1] whzm = 0.5*w[i,j ,k] + 0.5*w[i-1,j,k] elif axis == 1: uhxp = 0.5*u[i+1,j,k] + 0.5*u[i+1,j-1,k] uhxm = 0.5*u[i,j,k] + 0.5*u[i,j-1,k] vhyp = 0.5*v[i,j,k] + 0.5*v[i,j+1,k] vhym = 0.5*v[i,j,k] + 0.5*v[i,j-1,k] whzp = 0.5*w[i,j,k+1] + 0.5*w[i,j-1,k+1] whzm = 0.5*w[i,j ,k] + 0.5*w[i,j-1,k] elif axis == 2: uhxp = 0.5*u[i+1,j,k] + 0.5*u[i+1,j,k-1] uhxm = 0.5*u[i,j,k] + 0.5*u[i,j,k-1] vhyp = 0.5*v[i,j+1,k] + 0.5*v[i,j+1,k-1] vhym = 0.5*v[i,j,k] + 0.5*v[i,j,k-1] whzp = 0.5*w[i,j,k] + 0.5*w[i,j,k+1] whzm = 0.5*w[i,j,k] + 0.5*w[i,j,k-1] axp[c] = (max([-uhxp,0.0]) + nu/dx)*dy*dz axm[c] = (max([ uhxm,0.0]) + nu/dx)*dy*dz ayp[c] = (max([-vhyp,0.0]) + nu/dy)*dx*dz aym[c] = (max([ vhym,0.0]) + nu/dy)*dx*dz azp[c] = (max([-whzp,0.0]) + nu/dz)*dx*dy azm[c] = (max([ whzm,0.0]) + nu/dz)*dx*dy ac[c] = axp[c]+axm[c]+ayp[c]+aym[c]+azp[c]+azm[c] if axis == 0: b[c] = -dy*dz*( p[i,j,k] - p[i-1,j,k] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c] elif axis == 1: b[c] = -dx*dz*( p[i,j,k] - p[i,j-1,k] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c] elif axis == 2: b[c] = -dx*dy*( p[i,j,k] - p[i,j,k-1] ) + (1.0 - alphaU)/alphaU * ac[c] * V0[c] ac[c] /= alphaU act[c] = ac[c] if i == xmb : ac[c] += axm[c] elif i == xpb : ac[c] += axp[c] if j == ymb : ac[c] += aym[c] elif j == ypb : ac[c] += ayp[c] if k == zmb : ac[c] += azm[c] elif k == zpb : ac[c] += azp[c] b[c] += 2.0*azp[c]*Vbound def sor(u): u0 = np.copy(u) # deep copy for i in range(xs, nx): for j in range(ys, ny): for k in range(zs, nz): c = (i,j,k) xp,xm,yp,ym,zp,zm = (i+1,j,k),(i-1,j,k),(i,j+1,k),(i,j-1,k),(i,j,k+1),(i,j,k-1) if i == xmb : uxp ,uxm = u0[xp], zero elif i == xpb : uxp ,uxm = zero, u0[xm] else: uxp ,uxm = u0[xp], u0[xm] if j == ymb : uyp ,uym = u0[yp], zero elif j == ypb : uyp ,uym = zero , u0[ym] else: uyp ,uym = u0[yp], u0[ym] if k == zmb : uzp, uzm = u0[zp], zero elif k == zpb : uzp, uzm = zero, u0[zm] else: uzp ,uzm = u0[zp], u0[zm] ut = (axp[c]*uxp + axm[c]*uxm \ + ayp[c]*uyp + aym[c]*uym \ + azp[c]*uzp + azm[c]*uzm + b[c])/ac[c] u0[c] = ut return u0 unew = np.copy(Vdir) # deep copy for _ in range(100): utmp = sor(unew) err = np.abs(unew - utmp).max() unew = np.copy(utmp) if err < 10**(-3) : break # print(err) return unew, act
def caldp(p,u,v,w,acu,acv,acw,dx,dy,dz): zero = float(0.0) dp = np.zeros_like(p) ac = np.zeros_like(p) axp = np.zeros_like(p) axm = np.zeros_like(p) ayp = np.zeros_like(p) aym = np.zeros_like(p) azp = np.zeros_like(p) azm = np.zeros_like(p) b = np.zeros_like(p) nx, ny, nz = p.shape for i in range(nx): for j in range(ny): for k in range(nz): c = (i, j, k) axp[c] = zero if np.abs( acu[i+1,j, k ] ) < 10**(-8) else dy**2*dz**2/acu[i+1,j, k ] axm[c] = zero if np.abs( acu[i ,j, k ] ) < 10**(-8) else dy**2*dz**2/acu[i ,j, k ] ayp[c] = zero if np.abs( acv[i ,j+1,k ] ) < 10**(-8) else dx**2*dz**2/acv[i ,j+1,k ] aym[c] = zero if np.abs( acv[i ,j, k ] ) < 10**(-8) else dx**2*dz**2/acv[i ,j, k ] azp[c] = zero if np.abs( acw[i ,j, k+1] ) < 10**(-8) else dx**2*dy**2/acw[i ,j, k+1] azm[c] = zero if np.abs( acw[i ,j, k ] ) < 10**(-8) else dx**2*dy**2/acw[i ,j, k ] ac[c] = axp[c]+axm[c]+ayp[c]+aym[c]+azp[c]+azm[c] b[c] = - (u[i+1,j,k] - u[i,j,k])*dy*dz \ - (v[i,j+1,k] - v[i,j,k])*dx*dz \ - (w[i,j,k+1] - w[i,j,k])*dx*dy def sor(dp0, omega): p0 = np.copy(dp0) # deep copy for i in range(nx): for j in range(ny): for k in range(nz): c = (i,j,k) xp,xm,yp,ym,zp,zm = (i+1,j,k),(i-1,j,k),(i,j+1,k),(i,j-1,k),(i,j,k+1),(i,j,k-1) if i == 0: pxp ,pxm = p0[xp], zero elif i == nx-1: pxp ,pxm = zero, p0[xm] else: pxp ,pxm = p0[xp], p0[xm] if j == 0: pyp ,pym = p0[yp], zero elif j == ny-1: pyp ,pym = zero , p0[ym] else: pyp ,pym = p0[yp], p0[ym] if k == 0: pzp, pzm = p0[zp], zero elif k == nz-1: pzp, pzm = zero, p0[zm] else: pzp ,pzm = p0[zp], p0[zm] pt = (axp[c]*pxp + axm[c]*pxm \ + ayp[c]*pyp + aym[c]*pym \ + azp[c]*pzp + azm[c]*pzm + b[c])/ac[c] p0[c] += omega*(pt - p0[c]) # p0[c] = ut return p0 dpnew = np.copy(dp) # deep copy for nn in range(1000): dptmp = sor(dpnew, float(1.8)) err = np.abs(dpnew - dptmp).max() dpnew = np.copy(dptmp) if err < 10**(-6) : print('SOR N =', nn) break return dpnew
def SIMPLE(p,u,v,w,Uslip,Vslip,alphaU,alphaP,dx,dy,dz,nu): zero = float(0.0) nx, ny, nz = p.shape for nn in range(1000): u, acu = calupre(p,u,v,w,u,Uslip,alphaU,dx,dy,dz,nu,axis=0) v, acv = calupre(p,u,v,w,v,Vslip,alphaU,dx,dy,dz,nu,axis=1) w, acw = calupre(p,u,v,w,w,zero,alphaU,dx,dy,dz ,nu,axis=2) dp = caldp(p,u,v,w,acu,acv,acw,dx,dy,dz) p += alphaP*dp unew = np.copy(u) # deep copy vnew = np.copy(v) # deep copy wnew = np.copy(w) # deep copy for i in range(nx): for j in range(ny): for k in range(nz): if i !=0 : unew[i,j,k] = u[i,j,k] + (dp[i-1,j,k] - dp[i,j,k])*dy*dz/acu[i,j,k] if j !=0 : vnew[i,j,k] = v[i,j,k] + (dp[i,j-1,k] - dp[i,j,k])*dx*dz/acv[i,j,k] if k !=0 : wnew[i,j,k] = w[i,j,k] + (dp[i,j,k-1] - dp[i,j,k])*dx*dy/acw[i,j,k] err = max([np.abs(unew - u).max(), np.abs(vnew - v).max(), np.abs(wnew - w).max()]) print('N', nn, 'error =', err) u = np.copy(unew) # deep copy v = np.copy(vnew) # deep copy w = np.copy(wnew) # deep copy if err < 10**(-6) : break return p,u,v,w
プログラミングは以下を参考にしました。
二次元の図はHoloviewsで書いてます。 ソースはこちら
computational-sediment-hyd.hatenablog.jp
三次元の図はParaViewで書きました。 XMLベースのvtk形式に関する覚書 - Qiita を参考にxmlベースのvtkを作成しました。 ソースはこちら
話題のParaView Glance(web版ParaView)でレンダリングしてみました。
以下のサイトを参考にParaViewで作った図をvtkjsでエクスポートしました。以下ではhtmlも作成していますが上手くいきませんでした。
描画はParaView Glanceにアクセスしてvtkjsファイル(サンプル)をロードするだけです。 手順をgif化しました。
関連記事
久しぶりにParaviewを使うおうとインストールしたら不具合が出たので。
元ネタIntel GPU driver update 26.20.100.6708 on windows breaks Qt GUI --> FIXED by UPGRADING intel drivers (#19364) · Issues · ParaView / ParaView · GitLab が英語しかなかったので日本語でまとめておきます。
データファイルを読み込み「Apply」ボタンをクリックするとこんな感じにtoolbarがクラッシュしてしまいます。
「インテル グラフィックスの設定」のバージョンが古いと発生するエラーのようです。 そのため、「インテル グラフィックスの設定」を最新版にバージョンアップすれば解決します。
バージョンによってインターフェイスやアップデート方法が異なります。「インテル グラフィックスの設定」を開いて手動でアップデートして下さい。
「インテル グラフィックスの設定」はタスクバー右端のアイコンおよびコントロールパネルから開くことができます。
私の場合はバージョン27.20.100.8280にアップデートして解決しました。ちなみにプロキシ環境だと上手くいきませんでした。
アップデート後はこんな感じに無事解決しました。
HHKBに慣れてしまうとノートPCを使ったときにCaps Lockを間違えてしまうので。
5分でcapslockキーをctrlキーに変更する(windows10) | ジョイタスネット
アップデートされてないしレジストリいじるから危ないけど。
- 過去記事よりサンプル
ミャンマー某河川の河口部:ETOPO1 + geoviews - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning
単純な折れ線グラフをHoloviewsで書いてみます。比較のため、matplotlib版も書いておきます。 コードはjupyter上で書いています。なお、グラフはGitHubからGitHub Pagesを表示してiframeで取り込んでいます。
import holoviews as hv hv.extension('bokeh') g = hv.Curve([1,2,4],label='line1').options(width=400,height=300) * hv.Curve([1,4,8],label='line2') g
%matplotlib inline import matplotlib.pyplot as plt fig, ax = plt.subplots(figsize=(4.0, 3.0)) ax.plot([1,2,4],label='line1') ax.plot([1,4,8],label='line2') ax.legend()
個人的には静的なグラフなんていらないと思いますが、 いざ論文とかにまとめようと思うと静的なベクターのグラフを必要だったりします。(近い将来には変わると思いますが)
前述のとおり、Holoviewsはmatplotlibのラッパーでもあるのでバックエンドを切り替えるだけで簡単に静的なベクターのグラフを作ることができます。
import holoviews as hv hv.extension('matplotlib') g = hv.Curve([1,2,4],label='line1').options(aspect=1.333, fig_size=150) * hv.Curve([1,4,8],label='line2') hv.save(g,'holoviews.svg') g
縦横比の設定方法が、bokehとmatplotlibで異なるのでちょっと面倒ですが、if文とか書いておけば良いと思います。