周期境界は配列の循環シフトを使えばソースコードがスッキリするかも - 趣味で計算流砂水理 Computational Sediment Hydraulics for Fun Learning

循環シフトとは

配列の最初と最後を循環してシフトするものです。 pythonのnumpy.rollで書くとこんな感じになります。

import numpy as np

x = np.arange(5)
x # array([0, 1, 2, 3, 4]) 

np.roll(x,2)
# array([4, 0, 1, 2, 3])

np.roll(x,-1)
# array([1, 2, 3, 4, 0])

いまいち何に使うのかわかりませんが(ビット演算とかでしょうか)、これを使えば周期境界の離散式がスッキリかけると思います。

ちなみにfortranだとcshiftです。

python

numpy.roll

https://note.nkmk.me/python-numpy-roll/

fortran

cshift

https://www.sci.hokudai.ac.jp/~inaz/doc/B/Fortran90/node19.html

循環シフトを用いた周期境界のコーディング

以下の単純な移流方程式を対象とする。

$\begin{align} \dfrac{\partial u }{\partial t } &= - c\dfrac{\partial u }{\partial x } \\ c &= \rm{const} \end{align}$

計算条件として、 $c=1$ および $u$ の初期値を次式で与える。

$\begin{align} u_0 = \mathrm{e}^{-(x-ct)^2} \end{align}$

また、 $dt=0.01,dx=0.1$ とする。

x = np.linspace(-5,5,101, endpoint=True )
u = np.exp(-x**2)

c = 1.0
dt = 0.01
dx = 0.1

周期境界で $u$ の時間発展を中央差分で解くコードは次のようになる。

ut = np.zeros_like(u)
for i in range(len(u)):
    if i == len(u)-1:
        ut[i] = u[i] - dt * c * 0.5 * (u[0] - u[i-1])/dx
    else:
        ut[i] = u[i] - dt * c * 0.5 * (u[i+1] - u[i-1])/dx

pythonの場合、-1のインデックスは配列の末尾を示すため、配列の最後の計算だけをif文で分岐させる必要がある。

または、スライシングを使って次のようになる。

ut = np.zeros_like(u)
ut[1:len(u)-1] = u[1:len(u)-1] - dt * c * 0.5 * (u[2:len(u)] - u[0:len(u)-2])/dx

i = 0 
ut[i] = u[i] - dt * c * 0.5 * (u[i+1] - u[i-1])/dx

i = len(u)-1
ut[i] = u[i] - dt * c * 0.5 * (u[0] - u[i-1])/dx

境界部分は上手くスライシングできないため、別に書く必要がある。

次に、このコードを循環シフトを使って書くと次のようになる。

ut = u - dt * c * 0.5*( np.roll(u,-1) - np.roll(u,1) )/dx

if文が不要で、ループがいらないためインデックスも不要である。

roll関数が若干わかりにくいので好みはわかれるがワンライナー好きとしてはこちらの方が良い。

計算例

先程の移流方程式を循環シフトを使って中央差分と後退差分(風上差分)で計算してみた。

10秒後の計算結果を示す。なお、周期境界のため真値は初期値と同じになる。

import numpy as np
import holoviews as hv

hv.extension('bokeh')

x = np.linspace(-5,5,101, endpoint=True )
uini = np.exp(-x**2)
c = 1.0
dt = 0.01
dx = 0.1

#中央差分
u = np.copy(uini)
for i in range(1000):
    ut = u - dt * c * 0.5*( np.roll(u,-1) - np.roll(u,1) )/dx
    u = np.copy(ut)
uc = np.copy(ut)

#風上差分
u = np.copy(uini)
for i in range(1000):
    ut = u - dt * c * (u - np.roll(u,1))/dx
    u = np.copy(ut)
ub = np.copy(ut)

g = hv.Curve((x,uini),label='initial & true') * hv.Curve((x,uc), label='central') * hv.Curve((x,ub), label='upwind')
g.options(width=500, legend_position='right')

f:id:SedimentHydraulics:20201006193508p:plain