PMLのポーティング(20160315)

fdtd_3d_pml_v1_Cporting.py

#! coding:utf-8
"""
fdtd_3d.py
日本音響学会　サイエンスシリーズ14
「FDTD法で視る音の世界」　付録DVD
fdtd_3d_c_sjis.cのpython翻訳
所管
x->i
y->j
z->k
"""
import os
import numpy as np

# 変数の宣言
xmax = 5.000e0  # x軸解析領域 [m]
ymax = 5.000e0  # y軸解析領域 [m]
zmax = 5.000e0  # z軸解析領域 [m]
tmax = 2.000e-2  # 解析時間 [s]
dh = 5.000e-1  # 空間離散化幅 [m]
dt = 8.400e-5  # 時間離散化幅 [s]
c0 = 3.435e2  # 空気の音速 [m/s]
row0 = 1.205e0  # 空気の密度 [kg/m^3]
xdr = 2.000e0  # x軸音源位置 [m]
ydr = 3.000e0  # y軸音源位置 [m]
zdr = 2.500e0  # z軸音源位置 [m]
xon = 2.500e0  # 直方体x座標最小値 [m]
xox = 3.500e0  # 直方体x座標最大値 [m]
yon = 1.500e0  # 直方体y座標最小値 [m]
yox = 3.000e0  # 直方体y座標最大値 [m]
zon = 1.500e0  # 直方体z座標最小値 [m]
zox = 3.500e0  # 直方体z座標最大値 [m]
alpn = 0.200e0  # 直方体表面吸音率 [-]
m = 1.000e0  # ガウシアンパルス最大値 [m^3/s]
a = 2.000e6  # ガウシアンパルス係数 [-]
t0 = 3.000e-3  # ガウシアンパルス中心時間 [s]

pl = 16  # PML層数 [-]
pm = 4  # PML減衰係数テーパー乗数 [-]
emax = 1.200e0  # PML減衰係数最大値

fn = "out_3d"  # 出力ファイルネーム
df = 5  # 出力ファイルスキップ数

# Buffer
ex = np.zeros(pl + 1)
pmla = np.zeros(pl + 1)
pmlb = np.zeros(pl + 1)
pmlc = np.zeros(pl + 1)

# 諸定数の算出

# 解析範囲
ix = int(xmax / dh) + pl * 2
jx = int(ymax / dh) + pl * 2
kx = int(zmax / dh) + pl * 2
tx = int(tmax / dt)

# 直方体位置
ion = int(xon / dh) + pl
iox = int(xox / dh) + pl
jon = int(yon / dh) + pl
jox = int(yox / dh) + pl
kon = int(zon / dh) + pl
kox = int(zox / dh) + pl

# 加振点位置
idr = int(xdr / dh) + pl
jdr = int(ydr / dh) + pl
kdr = int(zdr / dh) + pl

# 加振時間
tdr = int((2.0 * t0) / dt)

# 体積弾性率
kp0 = row0 * c0 * c0

# 特性インピーダンス
z0 = row0 * c0

# 表面インピーダンス
if alpn != 0.0:
    zn = row0 * c0 * (1.0 + np.sqrt(1.0 - alpn)) / (1.0 - np.sqrt(1.0 - alpn))

# Courant数
clf = c0 * dt / dh

# 粒子速度用更新係数
vc = clf / z0

# 音圧用更新係数
pc = clf * z0

# PML用更新係数
for i in range(pl):
    i += 1
    ex[i] = emax * np.power(float(pl - i + 1) / float(pl), float(pm))

for i in range(pl):
    i += 1
    pmla[i] = (1.0 - ex[i]) / (1.0 + ex[i])
    pmlb[i] = clf / z0 / (1.0 + ex[i])
    pmlc[i] = clf * z0 / (1.0 + ex[i])

# メモリ格納
p = np.zeros((ix + 1, jx + 1, kx + 1))
px = np.zeros((ix + 1, jx + 1, kx + 1))
py = np.zeros((ix + 1, jx + 1, kx + 1))
pz = np.zeros((ix + 1, jx + 1, kx + 1))

vx = np.zeros((ix + 1, jx + 1, kx + 1))
vy = np.zeros((ix + 1, jx + 1, kx + 1))
vz = np.zeros((ix + 1, jx + 1, kx + 1))

q = np.zeros(tdr + 1)

print "p.shape", p.shape
print "px.shape", px.shape

# 音源波形の生成
for t in range(tdr):
    # t={1,2,...,tdr}
    t += 1
    q[t] = m * np.exp(-a * pow(float(t * dt - t0), 2.0))

# 時間ループ
tcount = 1
fcount = 0
txstep = float(tx) / 100.

print("Time Loop Start" + os.linesep)

for t in range(tx):
    # t = {1,2,...,tx}
    t += 1
    print t, tx

    # ----------------------------
    # 粒子速度(vx)の更新
    # ----------------------------

    # 左側のPML
    for i in range(pl):
        i += 1
        for j in range(jx):
            j += 1
            for k in range(kx):
                k += 1
                vx[i, j, k] = pmla[i] * vx[i, j, k] - pmlb[i] * (p[i + 1, j, k] - p[i, j, k])
                pass

    # 音響領域
    for i in range(pl + 1, ix - pl - 1 + 1):
        # print i, ix-pl-1+1
        for j in range(1, jx + 1):
            # print j, jx+1
            for k in range(1, kx + 1):
                vx[i, j, k] = vx[i, j, k] - vc * (p[i + 1, j, k] - p[i, j, k])

    # 左側PML
    for i in range(ix - pl, ix - 1 + 1):
        # for(i=ix-pl; i<=ix-1;i++)
        for j in range(1, jx + 1):
            for k in range(1, kx + 1):
                vx[i, j, k] = pmla[ix - i] * vx[i, j, k] - pmlb[ix - i] * (p[i + 1, j, k] - p[i, j, k])

    # -- (上) PMLの範囲(配列番号)がわかれば3行
    for i in range(1, ix + 1):
        for j in range(1, pl + 1):
            for k in range(1, kx + 1):
                vy[i, j, k] = pmla[j] * vy[i, j, k] - pmlb[j] * (p[i, j + 1, k] - p[i, j, k])

    for i in range(ix + 1):
        for j in range(pl + 1, jx - pl - 1 + 1):
            for k in range(1, kx + 1):
                vy[i, j, k] = vy[i, j, k] - vc * (p[i, j + 1, k] - p[i, j, k])

    for i in range(1, ix + 1):
        for j in range(jx - pl, jx - 1 + 1):
            for k in range(1, kx + 1):
                vy[i, j, k] = pmla[jx - j] * vy[i, j, k] - pmlb[jx - j] * (p[i, j + 1, k] - p[i, j, k])

    # 粒子速度(vz)の更新
    for i in range(1, ix + 1):
        for j in range(1, jx + 1):
            for k in range(1, pl + 1):
                vz[i, j, k] = pmla[k] * vz[i, j, k] - pmlb[k] * (p[i, j, k + 1] - p[i, j, k])

    for i in range(1, ix + 1):
        for j in range(1, jx + 1):
            for k in range(pl + 1, kx - pl - 1 + 1):
                vz[i, j, k] = vz[i, j, k] - vc * (p[i, j, k + 1] - p[i, j, k])

    for i in range(1, ix + 1):
        for j in range(1, jx + 1):
            for k in range(kx - pl, kx - 1 + 1):
                vz[i, j, k] = pmla[kx - k] * vz[i, j, k] - pmlb[kx - k] * (p[i, j, k + 1] - p[i, j, k])

    # 境界条件(vx)の計算
    for j in range(1, jx + 1):
        for k in range(1, kx + 1):
            vx[0, j, k] = 0.0
            vx[ix, j, k] = 0.0

    for j in range(jon, jox + 1):
        for k in range(kon, kox + 1):
            if alpn != 0.0:
                vx[ion - 1, j, k] = p[ion - 1, j, k] / zn
                vx[ion, j, k] = -p[iox + 1, j, k] / zn
            else:
                vx[ion - 1, j, k] = 0.0
                vx[iox, j, k] = 0.0

    # 境界条件(vy)の計算
    for k in range(1, kx + 1):
        for i in range(1, ix + 1):
            vy[i, 0, k] = 0.0
            vy[i, jx, k] = 0.0

    for k in range(kon, kox + 1):
        for i in range(ion, iox + 1):
            if alpn != 0.0:
                vy[i, jon - 1, k] = p[i, jon - 1, k] / zn
                vy[i, jox, k] = p[i, jox + 1, k] / zn
            else:
                vy[i, jon - 1, k] = 0.0
                vy[i, jox, k] = 0.0
    # 境界条件(vz)の計算
    for i in range(1, ix + 1):
        for j in range(1, jx + 1):
            vz[i, j, 0] = 0.0
            vz[i, j, kx] = 0.0
    for i in range(ion, iox + 1):
        for j in range(jon, jox + 1):
            if alpn != 0.0:
                vz[i, j, kon - 1] = p[i, j, kon - 1] / zn
                vz[i, j, kox] = -p[i, j, kox + 1] / zn
            else:
                vz[i, j, kon - 1] = 0.0
                vz[i, j, kox] = 0.0

    # 音圧(px)の更新
    for i in range(1, pl + 1):
        for j in range(1, jx + 1):
            for k in range(1, kx + 1):
                px[i, j, k] = pmla[i] * px[i, j, k] - pmlc[i] * (vx[i, j, k] - vx[i - 1, j, k])

    for i in range(pl + 1, ix - pl + 1):
        for j in range(1, jx + 1):
            for k in range(1, kx + 1):
                px[i, j, k] = px[i, j, k] - pc * (vx[i, j, k] - vx[i - 1, j, k])
                if (i == idr) and (j == jdr) and (k == kdr) and (t <= tdr):
                    px[i, j, k] = px[i, j, k] + dt * kp0 * q[t] / 3.0 / (dh * dh * dh)

    for i in range(ix - pl + 1, ix + 1):
        for j in range(1, jx + 1):
            for k in range(1, kx + 1):
                px[i, j, k] = pmla[ix - i + 1] * px[i, j, k] - pmlc[ix - i + 1] * (vx[i, j, k] - vx[i - 1, j, k])

    # 音圧(py)の更新
    for i in range(1, ix + 1):
        for j in range(1, pl + 1):
            for k in range(1, kx + 1):
                py[i, j, k] = pmla[j] * py[i, j, k] - pmlc[j] * (vy[i, j, k] - vy[i, j - 1, k])
    for i in range(1, ix + 1):
        for j in range(pl * 1, jx - pl + 1):
            for k in range(1, kx + 1):
                py[i, j, k] = p[i, j, k] - pc * (vy[i, j, k] - vy[i, j - 1, k])
                if (i == idr) and (j == jdr) and (k == kdr) and (t <= tdr):
                    py[i, j, k] = py[i, j, k] + dt * kp0 * q[t] / 3.0 / (dh * dh * dh)
    for i in range(1, ix + 1):
        for j in range(jx - pl + 1, jx + 1):
            for k in range(1, kx + 1):
                py[i, j, k] = pmla[jx - j + 1] * py[i, j, k] - pmlc[jx - j + 1] * (vy[i, j, k] - vy[i, j - 1, k])

    # 音圧(pz)の更新
    for i in range(1, ix + 1):
        for j in range(1, jx + 1):
            for k in range(1, pl + 1):
                pz[i, j, k] = pmla[k] * pz[i, j, k] - pmlc[k] * (vz[i, j, k] - vz[i, j, k - 1])

    for i in range(1, ix + 1):
        for j in range(1, jx + 1):
            for k in range(pl + 1, kx - pl + 1):
                pz[i, j, k] = pz[i, j, k] - pc * (vz[i, j, k] - vz[i, j, k - 1])

    for i in range(1, ix + 1):
        for j in range(1, jx + 1):
            for k in range(kx - pl + 1, kx + 1):
                pz[i, j, k] = pmla[kx - k + 1] * pz[i, j, k] - pmlc[kx - k + 1] * (vz[i, j, k] - vz[i, j, k - 1])

    # 音圧の合成
    p = px + py + pz


    os.linesep="\n"
    # 結果の出力
    fcount += 1
    filename = "%s_%05d.vtk" % (fn, fcount)
    dirname = fn
    if not os.path.exists(dirname):
        os.mkdir(dirname)
    filename = os.path.join(dirname,filename)
    with open(filename, 'w') as f:
        f.write("# vtk DataFile Version 2.0" + os.linesep)
        f.write(filename + os.linesep)
        f.write("ASCII" + os.linesep)
        f.write("DATASET STRUCTURED_POINTS" + os.linesep)
        f.write("DIMENSIONS %5d%5d%5d%s" % (ix - 2 * pl, jx - 2 * pl, kx - 2 * pl, os.linesep))
        f.write("ORIGIN %7.3f%7.3f%7.3f%s" % (0.0, 0.0, 0.0, os.linesep))
        f.write("SPACING %7.3f%7.3f%7.3f%s" % (dh, dh, dh, os.linesep))
        f.write("%s%10d%s" % ("POINT_DATA ", (ix - 2 * pl) * (jx - 2 * pl) * (kx - 2 * pl), os.linesep))
        f.write("SCALARS SoundPressure float" + os.linesep)
        f.write("%s%s" % ("LOOKUP_TABLE default", os.linesep))

        for k in range(pl + 1, kx - pl + 1):
            for j in range(pl + 1, jx - pl + 1):
                for i in range(pl + 1, ix - pl + 1):
                    if (i >= ion and i <= iox and j >= jon and j <= jox and k >= kon and k <= kox):
                        f.write("%20.10e%s" % (-100., os.linesep))
                    else:
                        f.write("%20.10e%s" % (p[i, j, k], os.linesep))

        f.write("%s%s" % ("VECTORS ParticleVelocity float", os.linesep))
        for k in range(pl + 1, kx - pl + 1):
            for j in range(pl + 1, jx - pl + 1):
                for i in range(pl + 1, ix - pl + 1):
                    if (i >= ion and i <= iox and j >= jon and j <= jox and k >= kon and k <= kox):
                        f.write("%20.10e%20.10e%20.10e%s" % (0., 0., 0., os.linesep))
                    else:

                        tmpx = (vx[i - 1, j, k] + vx[i, j, k]) / 2.0
                        tmpy = (vy[i, j - 1, k] + vy[i, j, k]) / 2.0
                        tmpz = (vz[i, j, k - 1] + vz[i, j, k]) / 2.0
                        f.write("%20.10e%20.10e%20.10e%s" % (tmpx, tmpy, tmpz, os.linesep))

    if (t >= int(float(tcount) * txstep)):
        print("%s%7.2f%s%s" % ("Completed", float(t) / float(tx) * 100., "%", os.linesep))
        tcount += 1

fdtd_3d_v2.py

#!coding:utf-8
"""
fdtd_3d_v1.py

FDTD
境界条件：インピーダンス境界
次元：3D

"""
import os
from time import sleep
from scipy import *
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import matplotlib.pyplot as plt

# ************************************** #
# * Paramater
# ************************************** #
X = 50
Y = 50
Z = 50

dx = 0.1
dy = 0.1
dz = 0.1
dt = 0.0001

Ro = 1.21
C = 343
K = Ro * C * C

# ************************************** #
# * Source
# ************************************** #
ang = arange(-pi, pi, 2 * pi / 50)
sig = cos(ang)
sig += 1

# 音響インピーダンス
Z0 = Ro * C

# ************************************** #
# * Buffer
# ************************************** #
P1 = zeros((X, Y, Z), "float64")
P2 = zeros((X, Y, Z), "float64")

Ux1 = zeros((X + 1, Y, Z), "float64")
Ux2 = zeros((X + 1, Y, Z), "float64")

Uy1 = zeros((X, Y + 1, Z), "float64")
Uy2 = zeros((X, Y + 1, Z), "float64")

Uz1 = zeros((X, Y, Z + 1), "float64")
Uz2 = zeros((X, Y, Z + 1), "float64")

Rec = []

# ************************************** #
# * Field Plot
# ************************************** #
plot_z = Z / 2

x = arange(0, dx * X, dx)
y = arange(0, dy * Y, dy)

xx, yy = meshgrid(y, x)

fig = plt.figure()
ax = axes3d.Axes3D(fig)

surf = ax.plot_surface(xx, yy, P2[:, :, plot_z], rstride=1, cstride=1, cmap=cm.jet)
oldsurf1 = surf
surf = ax.plot_surface(xx, yy, P2[:, :, plot_z / 2] + 1, rstride=1, cstride=1, cmap=cm.jet)
oldsurf2 = surf
fig.show()

# ************************************** #
# * Iteration
# ************************************** #
for t in range(220):
    print t
    if t < len(sig):
        P1[X / 2, Y / 2, Z / 2] += sig[t]

    # ** 運動方程式 ** #
    # Inner Elements
    Ux2[1:X, :, :] = Ux1[1:X, :, :] - dt / Ro / dx * (P1[1:X, :, :] - P1[:X - 1, :, :])
    Uy2[:, 1:Y, :] = Uy1[:, 1:Y, :] - dt / Ro / dy * (P1[:, 1:Y, :] - P1[:, :Y - 1, :])
    Uz2[:, :, 1:Z] = Uz1[:, :, 1:Z] - dt / Ro / dz * (P1[:, :, 1:Z] - P1[:, :, :Z - 1])

    # 振動解析により，構造振動ノードの変位から速度計算
    # write code ...

    # BC
    Ux2[0, :, :] = P1[0, :, :] / -Z0
    Ux2[-1, :, :] = P1[-1, :, :] / Z0

    Uy2[:, 0, :] = P1[:, 0, :] / -Z0
    Uy2[:, -1, :] = P1[:, -1, :] / Z0

    Uz2[:, :, 0] = P1[:, :, 0] / -Z0
    Uz2[:, :, -1] = P1[:, :, -1] / Z0


    # ** 連続の式 ** #
    P2[:X, :Y, :Z] = P1[:X, :Y, :Z] \
                     - K * dt / dx * (Ux2[1:X + 1, :, :] - Ux2[:X, :, :]) \
                     - K * dt / dy * (Uy2[:, 1:Y + 1, :] - Uy2[:, :Y, :]) \
                     - K * dt / dz * (Uz2[:, :, 1:Z + 1] - Uz2[:, :, :Z])

    # グラフプロット
    surf1 = ax.plot_surface(xx, yy, P2[:, :, plot_z], rstride=1, cstride=1, cmap=cm.jet)
    surf2 = ax.plot_surface(xx, yy, P2[:, :, plot_z / 2] + 1, rstride=1, cstride=1, cmap=cm.jet)
    ax.set_zlim3d(-1, 1)
    ax.collections.remove(oldsurf1)
    ax.collections.remove(oldsurf2)
    oldsurf1 = surf1
    oldsurf2 = surf2
    plt.pause(.01)

    # 変数の置き換え
    P1, P2 = P2, P1
    Ux1, Ux2 = Ux2, Ux1
    Uy1, Uy2 = Uy2, Uy1
    Uz1, Uz2 = Uz2, Uz1

    Rec.append(P2[X / 4, Y / 4, Z / 4])


    # output
    # 結果の出力

    os.linesep="\n"
    fcount = t
    fn = "out"
    filename = "out_%05d.vtk" % (fcount)
    dirname = os.path.basename(__file__)
    dirname, ext = os.path.splitext(dirname)
    dirname = "out_"+dirname
    if not os.path.exists(dirname):
        os.mkdir(dirname)
    filename = os.path.join(dirname,filename)

    # アウトプット
    with open(filename, 'w') as f:
        f.write("# vtk DataFile Version 2.0" + os.linesep)
        f.write(filename + os.linesep)
        f.write("ASCII" + os.linesep)
        f.write("DATASET STRUCTURED_POINTS" + os.linesep)
        f.write("DIMENSIONS %5d%5d%5d%s" % (X, Y, Z, os.linesep))
        f.write("ORIGIN %7.3f%7.3f%7.3f%s" % (0.0, 0.0, 0.0, os.linesep))
        f.write("SPACING %7.3f%7.3f%7.3f%s" % (dx, dy, dz, os.linesep))
        f.write("%s%10d%s" % ("POINT_DATA ", X*Y*Z, os.linesep))
        f.write("SCALARS SoundPressure float" + os.linesep)
        f.write("%s%s" % ("LOOKUP_TABLE default", os.linesep))

        for k in range(Z):
            for j in range(Y):
                for i in range(X):
                    f.write("%20.10e%s" % (P2[i, j, k], os.linesep))


        f.write("%s%s" % ("VECTORS ParticleVelocity float", os.linesep))
        for k in range(Z):
            for j in range(Y):
                for i in range(X):
                        tmpx = (Ux2[i - 1, j, k] + Ux2[i, j, k]) / 2.0
                        tmpy = (Uy2[i, j - 1, k] + Uy2[i, j, k]) / 2.0
                        tmpz = (Uz2[i, j, k - 1] + Uz2[i, j, k]) / 2.0
                        f.write("%20.10e%20.10e%20.10e%s" % (tmpx, tmpy, tmpz, os.linesep))


fig = plt.figure()

plt.plot(Rec / max(Rec))
plt.plot(sig / max(sig))
plt.show()