fdtd_3d_v1.py

fdtd_3d_v1.py

#!coding:utf-8
"""
fdtd_3d_v1.py
OpenAcousticsのサンプルコード

"""
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 = 10
Y = 10
Z = 10

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

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


    # output
    filename = "res_000.pp"
    basename = os.path.splitext(filename)[0]

    with open(basename + ".vtk", "w") as file:

        file.write("# vtk DataFile Version 2.0\n")
        file.write(basename + "\n")
        file.write("ASCII\n")
        file.write("DATASET UNSTRUCTURED_GRID\n")

        file.write("POINTS %d float\n" % len(x))





fig = plt.figure()
plt.plot(Rec / max(Rec))
plt.plot(sig / max(sig))
plt.show()