Skip to content

Instantly share code, notes, and snippets.

@peace098beat

peace098beat/fdtd_2d_c_sjis.c

Last active February 8, 2020 13:24
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save peace098beat/6cb8d051919d19511d7c to your computer and use it in GitHub Desktop.
Save peace098beat/6cb8d051919d19511d7c to your computer and use it in GitHub Desktop.
Download ZIP
[FDTD] FDTDによる音の可視化 C->pythonにポーティング
Raw
fdtd_2d_c_sjis.c
// 二次元音響FDTDサンプル for C (by Masahiro TOYODA)
// gnu gcc -> gcc fdtd_2d_c_sjis.c -lm -O3 -fopenmp -o fdtd_2d_c_sjis.exe
// intel c -> icc fdtd_2d_c_sjis.c /O3 /Qopenmp
//■■■ヘッダーの読み込み■■■
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
//■■■定数の宣言■■■
#define xmax 5.000e0 // x軸解析領域 [m]
#define ymax 5.000e0 // y軸解析領域 [m]
#define tmax 2.000e-2 // 解析時間 [s]
#define dh 5.000e-2 // 空間離散化幅 [m]
#define dt 1.000e-4 // 時間離散化幅幅 [s]
#define c0 3.435e2 // 空気の音速 [m/s]
#define row0 1.205e0 // 空気の密度 [kg/m^3]
#define xdr 2.000e0 // x軸音源位置 [m]
#define ydr 3.000e0 // y軸音源位置 [m]
#define xon 2.500e0 // 直方体x座標最小値 [m]
#define xox 3.500e0 // 直方体x座標最大値 [m]
#define yon 1.500e0 // 直方体y座標最小値 [m]
#define yox 3.000e0 // 直方体y座標最大値 [m]
#define alpn 0.200e0 // 直方体表面吸音率 [-]
#define m 1.000e0 // ガウシアンパルス最大値 [m^3/s]
#define a 2.000e6 // ガウシアンパルス係数 [-]
#define t0 3.000e-3 // ガウシアンパルス中心時間 [s]
#define pl 16 // PML層数 [-]
#define pm 4 // PML減衰係数テーパー乗数 [-]
#define emax 1.500e0 // PML減衰係数最大値
#define fn "test_2d" // 出力ファイルネーム
#define df 5 // 出力ファイルスキップ数
//■■■変数の宣言■■■
double **p;
double **px;
double **py;
double **vx;
double **vy;
double *q;
double ex[pl+1];
double pmla[pl+1];
double pmlb[pl+1];
double pmlc[pl+1];
int t, i, j, ix, jx, ion, iox, jon, jox, tx, idr, jdr, tdr, tcount, fcount;
double kp0, z0, zn, clf, ex0, pc, vc, txstep;
char str[256];
FILE *fp;
int main(void)
{
//■■■諸定数の算出■■■
//■解析範囲■
ix = (int)(xmax / dh) + pl * 2;
jx = (int)(ymax / 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;
//■加振点位置■
idr = (int)(xdr / dh) + pl;
jdr = (int)(ydr / dh) + pl;
//■加振時間■
tdr = (int)((2.0 * t0) / dt);
//■体積弾性率■
kp0 = row0 * c0 * c0;
//■特性インピーダンス■
z0 = row0 * c0;
//■表面インピーダンス■
if(alpn != 0.0)
{
zn = row0 * c0 * (1.0 + sqrt(1.0 - alpn)) / (1.0 - sqrt(1.0 - alpn));
}
//■Courant数■
clf = c0 * dt / dh;
//■粒子速度用更新係数■
vc = clf / z0;
//■音圧用更新係数■
pc = clf * z0;
//■PML用更新係数■
#pragma omp parallel for private(i)
for(i = 1; i <= pl; i++)
{
ex[i] = emax * pow((double)(pl - i + 1) / (double)pl, (double)pm);
}
#pragma omp parallel for private(i)
for(i = 1; i <= pl; i++)
{
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 = malloc(sizeof(double*) * (ix+1));
for(i = 0; i <= ix; i++)
{
p[i] = malloc(sizeof(double) * (jx+1));
}
px = malloc(sizeof(double*) * (ix+1));
for(i = 0; i <= ix; i++)
{
px[i] = malloc(sizeof(double) * (jx+1));
}
py = malloc(sizeof(double*) * (ix+1));
for(i = 0; i <= ix; i++)
{
py[i] = malloc(sizeof(double) * (jx+1));
}
vx = malloc(sizeof(double*) * (ix+1));
for(i = 0; i <= ix; i++)
{
vx[i] = malloc(sizeof(double) * (jx+1));
}
vy = malloc(sizeof(double*) * (ix+1));
for(i = 0; i <= ix; i++)
{
vy[i] = malloc(sizeof(double) * (jx+1));
}
q = malloc(sizeof(double) * (tdr+1));
//■■■変数の初期化■■■
#pragma omp parallel for private(i, j)
for(i = 0; i <= ix; i++)
{
for(j = 0; j <= jx; j++)
{
p[i][j] = 0.0;
px[i][j] = 0.0;
py[i][j] = 0.0;
vx[i][j] = 0.0;
vy[i][j] = 0.0;
}
}
#pragma omp parallel for private(t)
for(t = 1; t <= tdr; t++)
{
q[t] = 0.0;
}
//■■■音源波形の作成■■■
#pragma omp parallel for private(t)
for(t = 1; t <= tdr; t++)
{
q[t] = m * exp(-a * pow((double)t * dt - t0, 2.0));
}
//■■■時間ループ■■■
tcount = 1;
fcount = 0;
txstep = (double)tx / 100.0;
printf("%s\n", "Time Loop Start");
for(t = 1; t <= tx; t++)
{
//■粒子速度(vx)の更新■
#pragma omp parallel for private(i, j)
for(i = 1; i <= pl; i++)
{
for(j = 1; j <= jx; j++)
{
vx[i][j] = pmla[i] * vx[i][j] - pmlb[i] * (p[i+1][j] - p[i][j]);
}
}
#pragma omp parallel for private(i, j)
for(i = pl + 1; i <= ix - pl - 1; i++)
{
for(j = 1; j <= jx; j++)
{
vx[i][j] = vx[i][j] - vc * (p[i+1][j] - p[i][j]);
}
}
#pragma omp parallel for private(i, j)
for(i = ix - pl; i <= ix - 1; i++)
{
for(j = 1; j <= jx; j++)
{
vx[i][j] = pmla[ix-i] * vx[i][j] - pmlb[ix-i] * (p[i+1][j] - p[i][j]);
}
}
//■粒子速度(vy)の更新■
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= pl; j++)
{
vy[i][j] = pmla[j] * vy[i][j] - pmlb[j] * (p[i][j+1] - p[i][j]);
}
}
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = pl + 1; j <= jx - pl - 1; j++)
{
vy[i][j] = vy[i][j] - vc * (p[i][j+1] - p[i][j]);
}
}
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = jx - pl; j <= jx - 1; j++)
{
vy[i][j] = pmla[jx-j] * vy[i][j] - pmlb[jx-j] * (p[i][j+1] - p[i][j]);
}
}
//■境界条件(vx)の計算■
#pragma omp parallel for private(j)
for(j = 1; j <= jx; j++)
{
vx[0][j] = 0.0;
vx[ix][j] = 0.0;
}
#pragma omp parallel for private(j)
for(j = jon; j <= jox; j++)
{
if(alpn != 0.0)
{
vx[ion-1][j] = p[ion-1][j] / zn;
vx[iox][j] = -p[iox+1][j] / zn;
}
else
{
vx[ion-1][j] = 0.0;
vx[iox][j] = 0.0;
}
}
//■境界条件(vy)の計算■
#pragma omp parallel for private(i)
for(i = 1; i <= ix; i++)
{
vy[i][0] = 0.0;
vy[i][jx] = 0.0;
}
#pragma omp parallel for private(i)
for(i = ion; i <= iox; i++)
{
if(alpn != 0.0)
{
vy[i][jon-1] = p[i][jon-1] / zn;
vy[i][jox] = -p[i][jox+1] / zn;
}
else
{
vy[i][jon-1] = 0.0;
vy[i][jox] = 0.0;
}
}
//■音圧(px)の更新■
#pragma omp parallel for private(i, j)
for(i = 1; i <= pl; i++)
{
for(j = 1; j <= jx; j++)
{
px[i][j] = pmla[i] * px[i][j] - pmlc[i] * (vx[i][j] - vx[i-1][j]);
}
}
#pragma omp parallel for private(i, j)
for(i = pl + 1; i <= ix - pl; i++)
{
for(j = 1; j <= jx; j++)
{
px[i][j] = px[i][j] - pc * (vx[i][j] - vx[i-1][j]);
if(i == idr && j == jdr && t <= tdr)
{
px[i][j] = px[i][j] + dt * kp0 * q[t] / 2.0 / (dh * dh);
}
}
}
#pragma omp parallel for private(i, j)
for(i = ix - pl + 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
px[i][j] = pmla[ix-i+1] * px[i][j] - pmlc[ix-i+1] * (vx[i][j] - vx[i-1][j]);
}
}
//■音圧(py)の更新■
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= pl; j++)
{
py[i][j] = pmla[j] * py[i][j] - pmlc[j] * (vy[i][j] - vy[i][j-1]);
}
}
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = pl + 1; j <= jx - pl; j++)
{
py[i][j] = py[i][j] - pc * (vy[i][j] - vy[i][j-1]);
if(i == idr && j == jdr && t <= tdr)
{
py[i][j] = py[i][j] + dt * kp0 * q[t] / 2.0 / (dh * dh);
}
}
}
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = jx - pl + 1; j <= jx; j++)
{
py[i][j] = pmla[jx-j+1] * py[i][j] - pmlc[jx-j+1] * (vy[i][j] - vy[i][j-1]);
}
}
//■音圧の合成■
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
p[i][j] = px[i][j] + py[i][j];
}
}
//■結果の出力■
if(t % df == 0)
{
fcount++;
sprintf(str, "%s_%05d.vtk", fn, fcount);
if((fp = fopen(str, "w")) != NULL)
{
fprintf(fp, "%s\n", "# vtk DataFile Version 2.0");
fprintf(fp, "%s\n", str);
fprintf(fp, "%s\n", "ASCII");
fprintf(fp, "%s\n", "DATASET STRUCTURED_POINTS");
fprintf(fp, "%s%5d%5d%5d\n", "DIMENSIONS ", ix - 2 * pl, jx - 2 * pl, 1);
fprintf(fp, "%s%7.3f%7.3f%7.3f\n", "ORIGIN ", 0.0, 0.0, 0.0);
fprintf(fp, "%s%7.3f%7.3f%7.3f\n", "SPACING ", dh, dh, 0.0);
fprintf(fp, "%s%10d\n", "POINT_DATA ", (ix - 2 * pl) * (jx - 2 * pl));
fprintf(fp, "%s\n", "SCALARS SoundPressure float");
fprintf(fp, "%s\n", "LOOKUP_TABLE default");
for(j = pl + 1; j <= jx - pl; j++)
{
for(i = pl + 1; i <= ix - pl; i++)
{
if(i >= ion && i <= iox && j >= jon && j <= jox)
{
fprintf(fp, "%20.10e\n", -100.0);
}
else
{
fprintf(fp, "%20.10e\n", p[i][j]);
}
}
}
fprintf(fp,"%s\n", "VECTORS ParticleVelocity float");
for(j = pl + 1; j <= jx - pl; j++)
{
for(i = pl + 1; i <= ix - pl; i++)
{
if(i >= ion && i <= iox && j >= jon && j <= jox)
{
fprintf(fp, "%20.10e%20.10e%20.10e\n", 0.0, 0.0, 0.0);
}
else
{
fprintf(fp, "%20.10e%20.10e%20.10e\n", (vx[i-1][j] + vx[i][j]) / 2.0, (vy[i][j-1] + vy[i][j]) / 2.0, 0.0);
}
}
}
}
fclose(fp);
}
//■進捗の出力■
if(t >= (int)((double)tcount * txstep))
{
printf("%s%7.2f%s\n", "Completed ", (double)t / (double)tx * 100.0, "%");
tcount++;
}
}
//■■■メモリの破棄■■■
free(p);
free(px);
free(py);
free(vx);
free(vy);
free(q);
return 0;
}
Raw
fdtd_2d_matlab_sjis.m
% 最も簡単な二次元音響FDTDサンプル for MATLAB (by Yoshiki NAGATANI)
% モデルの各定数
nx = 300; % X 方向の空間セル数 [pixels]
ny = 200; % Y 方向の空間セル数 [pixels]
dx = 10.0e-3; % 空間刻み [m]
dt = 20.0e-6; % 時間刻み [s]
nmax = 1000; % 計算ステップ数 [回]
savestep = 10; % 保存間隔（ステップ数）
% 媒質の各定数
rho = 1.293; % 媒質の密度ρ [kg/m^3]
kappa = 142.0e3; % 媒質の体積弾性率κ [Pa]
% 点音源の各定数
sig_freq = 1.0e3; % 印可波形の周波数 [Hz]
sig_amp = 5000.0; % 印可波形の振幅
sig_x = 100; % 点音源の x 座標 [pixels] ( 0 < sig_x <= nx)
sig_y = 80; % 点音源の y 座標 [pixels] ( 0 < sig_y <= ny)
% 音場の初期化
Vx = zeros(nx+1,ny ); % x方向粒子速度 [m/s]
Vy = zeros(nx, ny+1); % y方向粒子速度 [m/s]
P = zeros(nx, ny ); % 音圧 [Pa]
% ウィンドウの準備
h = figure();
colormap(jet(64));
% メインループ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for n = 0:1:nmax,
% 粒子速度の更新
Vx(2:nx,:) = Vx(2:nx,:) - dt / (rho * dx) * ( P(2:nx,:) - P(1:nx-1,:) );
Vy(:,2:ny) = Vy(:,2:ny) - dt / (rho * dx) * ( P(:,2:ny) - P(:,1:ny-1) );
% 音圧の更新
P(1:nx,1:ny) = P(1:nx,1:ny) - ( kappa * dt / dx ) * ( ( Vx(2:nx+1,:) - Vx(1:nx,:) ) + ( Vy(:,2:ny+1) - Vy(:,1:ny) ) );
% 点音源（正弦波×１波 with 二乗余弦関数）
if n < (1.0/sig_freq)/dt,
sig = sig_amp * (1.0-cos(2.0*pi*sig_freq*n*dt))/2.0 * sin(2.0*pi*sig_freq*n*dt);
P(sig_x,sig_y) = P(sig_x,sig_y) + sig;
end
% 音場の可視化（音圧）
image(abs(P)'); % 行列の画像表示（絶対値を表示する場合）：↓いずれかのみ
% image((32+P)'); % 行列の画像表示（正負の値を表示する場合）：↑いずれかのみ
axis equal; axis tight; % 軸の調整
title(sprintf('%7.3f ms [ %4d/%4d ] ([Ctrl]-[c] to STOP)', n*dt*1.0e3, n, nmax));
xlabel('x'); ylabel('y');
drawnow; % ウィンドウの表示内容を即座に更新する
% 画像ファイル保存（画像ファイルに書き出す場合は以下の４行をコメントアウト）
% if mod(n, savestep) == 0
% savename = sprintf('image_%04d.png',n);
% saveas(h,savename,'png');
% end
end
Raw
fdtd_3d_c_sjis.c
// 三次元音響FDTDサンプル for C (by Masahiro TOYODA)
// gnu gcc -> gcc fdtd_3d_c_sjis.c -lm -O3 -fopenmp -o fdtd_3d_c_sjis.exe
// intel c -> icc fdtd_3d_c_sjis.c /O3 /Qopenmp
//■■■ヘッダーの読み込み■■■
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
//■■■定数の宣言■■■
#define xmax 5.000e0 // x軸解析領域 [m]
#define ymax 5.000e0 // y軸解析領域 [m]
#define zmax 5.000e0 // z軸解析領域 [m]
#define tmax 2.000e-2 // 解析時間 [s]
#define dh 5.000e-2 // 空間離散化幅 [m]
#define dt 8.400e-5 // 時間離散化幅 [s]
#define c0 3.435e2 // 空気の音速 [m/s]
#define row0 1.205e0 // 空気の密度 [kg/m^3]
#define xdr 2.000e0 // x軸音源位置 [m]
#define ydr 3.000e0 // y軸音源位置 [m]
#define zdr 2.500e0 // z軸音源位置 [m]
#define xon 2.500e0 // 直方体x座標最小値 [m]
#define xox 3.500e0 // 直方体x座標最大値 [m]
#define yon 1.500e0 // 直方体y座標最小値 [m]
#define yox 3.000e0 // 直方体y座標最大値 [m]
#define zon 1.500e0 // 直方体z座標最小値 [m]
#define zox 3.500e0 // 直方体z座標最大値 [m]
#define alpn 0.200e0 // 直方体表面吸音率 [-]
#define m 1.000e0 // ガウシアンパルス最大値 [m^3/s]
#define a 2.000e6 // ガウシアンパルス係数 [-]
#define t0 3.000e-3 // ガウシアンパルス中心時間 [s]
#define pl 16 // PML層数 [-]
#define pm 4 // PML減衰係数テーパー乗数 [-]
#define emax 1.200e0 // PML減衰係数最大値
#define fn "test_3d" // 出力ファイルネーム
#define df 5 // 出力ファイルスキップ数
//■■■変数の宣言■■■
double ***p;
double ***px;
double ***py;
double ***pz;
double ***vx;
double ***vy;
double ***vz;
double *q;
double ex[pl+1];
double pmla[pl+1];
double pmlb[pl+1];
double pmlc[pl+1];
int t, i, j, k, ix, jx, kx, ion, iox, jon, jox, kon, kox, tx, idr, jdr, kdr, tdr, tcount, fcount;
double kp0, z0, zn, clf, ex0, pc, vc, txstep;
char str[256];
FILE *fp;
int main(void)
{
//■■■諸定数の算出■■■
//■解析範囲■
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 + sqrt(1.0 - alpn)) / (1.0 - sqrt(1.0 - alpn));
}
//■Courant数■
clf = c0 * dt / dh;
//■粒子速度用更新係数■
vc = clf / z0;
//■音圧用更新係数■
pc = clf * z0;
//■PML用更新係数■
#pragma omp parallel for private(i)
for(i = 1; i <= pl; i++)
{
ex[i] = emax * pow((double)(pl - i + 1) / (double)pl, (double)pm);
}
#pragma omp parallel for private(i)
for(i = 1; i <= pl; i++)
{
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 = malloc(sizeof(double**) * (ix+1));
for(i = 0; i <= ix; i++)
{
p[i] = malloc(sizeof(double*) * (jx+1));
for(j = 0; j <= jx; j++)
{
p[i][j] = malloc(sizeof(double) * (kx+1));
}
}
px = malloc(sizeof(double**) * (ix+1));
for(i = 0; i <= ix; i++)
{
px[i] = malloc(sizeof(double*) * (jx+1));
for(j = 0; j <= jx; j++)
{
px[i][j] = malloc(sizeof(double) * (kx+1));
}
}
py = malloc(sizeof(double**) * (ix+1));
for(i = 0; i <= ix; i++)
{
py[i] = malloc(sizeof(double*) * (jx+1));
for(j = 0; j <= jx; j++)
{
py[i][j] = malloc(sizeof(double) * (kx+1));
}
}
pz = malloc(sizeof(double**) * (ix+1));
for(i = 0; i <= ix; i++)
{
pz[i] = malloc(sizeof(double*) * (jx+1));
for(j = 0; j <= jx; j++)
{
pz[i][j] = malloc(sizeof(double) * (kx+1));
}
}
vx = malloc(sizeof(double**) * (ix+1));
for(i = 0; i <= ix; i++)
{
vx[i] = malloc(sizeof(double*) * (jx+1));
for(j = 0; j <= jx; j++)
{
vx[i][j] = malloc(sizeof(double) * (kx+1));
}
}
vy = malloc(sizeof(double**) * (ix+1));
for(i = 0; i <= ix; i++)
{
vy[i] = malloc(sizeof(double*) * (jx+1));
for(j = 0; j <= jx; j++)
{
vy[i][j] = malloc(sizeof(double) * (kx+1));
}
}
vz = malloc(sizeof(double**) * (ix+1));
for(i = 0; i <= ix; i++)
{
vz[i] = malloc(sizeof(double*) * (jx+1));
for(j = 0; j <= jx; j++)
{
vz[i][j] = malloc(sizeof(double) * (kx+1));
}
}
q = malloc(sizeof(double) * (tdr+1));
//■■■変数の初期化■■■
#pragma omp parallel for private(i, j, k)
for(i = 0; i <= ix; i++)
{
for(j = 0; j <= jx; j++)
{
for(k = 0; k <= kx; k++)
{
p[i][j][k] = 0.0;
px[i][j][k] = 0.0;
py[i][j][k] = 0.0;
pz[i][j][k] = 0.0;
vx[i][j][k] = 0.0;
vy[i][j][k] = 0.0;
vz[i][j][k] = 0.0;
}
}
}
#pragma omp parallel for private(t)
for(t = 1; t <= tdr; t++)
{
q[t] = 0.0;
}
//■■■音源波形の作成■■■
#pragma omp parallel for private(t)
for(t = 1; t <= tdr; t++)
{
q[t] = m * exp(-a * pow((double)t * dt - t0, 2.0));
}
//■■■時間ループ■■■
tcount = 1;
fcount = 0;
txstep = (double)tx / 100.0;
printf("%s\n", "Time Loop Start");
for(t = 1; t <= tx; t++)
{
//■粒子速度(vx)の更新■
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= pl; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
vx[i][j][k] = pmla[i] * vx[i][j][k] - pmlb[i] * (p[i+1][j][k] - p[i][j][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = pl + 1; i <= ix - pl - 1; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
vx[i][j][k] = vx[i][j][k] - vc * (p[i+1][j][k] - p[i][j][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = ix - pl; i <= ix - 1; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
vx[i][j][k] = pmla[ix-i] * vx[i][j][k] - pmlb[ix-i] * (p[i+1][j][k] - p[i][j][k]);
}
}
}
//■粒子速度(vy)の更新■
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= pl; j++)
{
for(k = 1; k <= kx; k++)
{
vy[i][j][k] = pmla[j] * vy[i][j][k] - pmlb[j] * (p[i][j+1][k] - p[i][j][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = pl + 1; j <= jx - pl - 1; j++)
{
for(k = 1; k <= kx; k++)
{
vy[i][j][k] = vy[i][j][k] - vc * (p[i][j+1][k] - p[i][j][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = jx - pl; j <= jx - 1; j++)
{
for(k = 1; k <= kx; k++)
{
vy[i][j][k] = pmla[jx-j] * vy[i][j][k] - pmlb[jx-j] * (p[i][j+1][k] - p[i][j][k]);
}
}
}
//■粒子速度(vz)の更新■
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= pl; k++)
{
vz[i][j][k] = pmla[k] * vz[i][j][k] - pmlb[k] * (p[i][j][k+1] - p[i][j][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = pl + 1; k <= kx - pl - 1; k++)
{
vz[i][j][k] = vz[i][j][k] - vc * (p[i][j][k+1] - p[i][j][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = kx - pl; k <= kx - 1; k++)
{
vz[i][j][k] = pmla[kx-k] * vz[i][j][k] - pmlb[kx-k] * (p[i][j][k+1] - p[i][j][k]);
}
}
}
//■境界条件(vx)の計算■
#pragma omp parallel for private(j, k)
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
vx[0][j][k] = 0.0;
vx[ix][j][k] = 0.0;
}
}
#pragma omp parallel for private(j, k)
for(j = jon; j <= jox; j++)
{
for(k = kon; k <= kox; k++)
{
if(alpn != 0.0)
{
vx[ion-1][j][k] = p[ion-1][j][k] / zn;
vx[iox][j][k] = -p[iox+1][j][k] / zn;
}
else
{
vx[ion-1][j][k] = 0.0;
vx[iox][j][k] = 0.0;
}
}
}
//■境界条件(vy)の計算■
#pragma omp parallel for private(k, i)
for(k = 1; k <= kx; k++)
{
for(i = 1; i <= ix; i++)
{
vy[i][0][k] = 0.0;
vy[i][jx][k] = 0.0;
}
}
#pragma omp parallel for private(k, i)
for(k = kon; k <= kox; k++)
{
for(i = ion; i <= iox; i++)
{
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)の計算■
#pragma omp parallel for private(i, j)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
vz[i][j][0] = 0.0;
vz[i][j][kx] = 0.0;
}
}
#pragma omp parallel for private(i, j)
for(i = ion; i <= iox; i++)
{
for(j = jon; j <= jox; j++)
{
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)の更新■
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= pl; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
px[i][j][k] = pmla[i] * px[i][j][k] - pmlc[i] * (vx[i][j][k] - vx[i-1][j][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = pl + 1; i <= ix - pl; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
px[i][j][k] = px[i][j][k] - pc * (vx[i][j][k] - vx[i-1][j][k]);
if((i == idr) && (j == jdr) && (k == kdr) && (t <= tdr))
{
px[i][j][k] = px[i][j][k] + dt * kp0 * q[t] / 3.0 / (dh * dh * dh);
}
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = ix - pl + 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
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)の更新■
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= pl; j++)
{
for(k = 1; k <= kx; k++)
{
py[i][j][k] = pmla[j] * py[i][j][k] - pmlc[j] * (vy[i][j][k] - vy[i][j-1][k]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = pl + 1; j <= jx - pl; j++)
{
for(k = 1; k <= kx; k++)
{
py[i][j][k] = py[i][j][k] - pc * (vy[i][j][k] - vy[i][j-1][k]);
if((i == idr) && (j == jdr) && (k == kdr) && (t <= tdr))
{
py[i][j][k] = py[i][j][k] + dt * kp0 * q[t] / 3.0 / (dh * dh * dh);
}
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = jx - pl + 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
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)の更新■
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= pl; k++)
{
pz[i][j][k] = pmla[k] * pz[i][j][k] - pmlc[k] * (vz[i][j][k] - vz[i][j][k-1]);
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = pl + 1; k <= kx - pl; k++)
{
pz[i][j][k] = pz[i][j][k] - pc * (vz[i][j][k] - vz[i][j][k-1]);
if((i == idr) && (j == jdr) && (k == kdr) && (t <= tdr))
{
pz[i][j][k] = pz[i][j][k] + dt * kp0 * q[t] / 3.0 / (dh * dh * dh);
}
}
}
}
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = kx - pl + 1; k <= kx; k++)
{
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]);
}
}
}
//■音圧の合成■
#pragma omp parallel for private(i, j, k)
for(i = 1; i <= ix; i++)
{
for(j = 1; j <= jx; j++)
{
for(k = 1; k <= kx; k++)
{
p[i][j][k] = px[i][j][k] + py[i][j][k] + pz[i][j][k];
}
}
}
//■結果の出力■
if(t % df == 0)
{
fcount++;
sprintf(str, "%s_%05d.vtk", fn, fcount);
if((fp = fopen(str, "w")) != NULL)
{
fprintf(fp, "%s\n", "# vtk DataFile Version 2.0");
fprintf(fp, "%s\n", str);
fprintf(fp, "%s\n", "ASCII");
fprintf(fp, "%s\n", "DATASET STRUCTURED_POINTS");
fprintf(fp, "%s%5d%5d%5d\n", "DIMENSIONS ", ix - 2 * pl, jx - 2 * pl, kx - 2 * pl);
fprintf(fp, "%s%7.3f%7.3f%7.3f\n", "ORIGIN ", 0.0, 0.0, 0.0);
fprintf(fp, "%s%7.3f%7.3f%7.3f\n", "SPACING ", dh, dh, dh);
fprintf(fp, "%s%10d\n", "POINT_DATA ", (ix - 2 * pl) * (jx - 2 * pl) * (kx - 2 * pl));
fprintf(fp, "%s\n", "SCALARS SoundPressure float");
fprintf(fp, "%s\n", "LOOKUP_TABLE default");
for(k = pl + 1; k <= kx - pl; k++)
{
for(j = pl + 1; j <= jx - pl; j++)
{
for(i = pl + 1; i <= ix - pl; i++)
{
if(i >= ion && i <= iox && j >= jon && j <= jox && k >= kon && k <= kox)
{
fprintf(fp, "%20.10e\n", -100.0);
}
else
{
fprintf(fp, "%20.10e\n", p[i][j][k]);
}
}
}
}
fprintf(fp,"%s\n", "VECTORS ParticleVelocity float");
for(k = pl + 1; k <= kx - pl; k++)
{
for(j = pl + 1; j <= jx - pl; j++)
{
for(i = pl + 1; i <= ix - pl; i++)
{
if(i >= ion && i <= iox && j >= jon && j <= jox && k >= kon && k <= kox)
{
fprintf(fp, "%20.10e%20.10e%20.10e\n", 0.0, 0.0, 0.0);
}
else
{
fprintf(fp, "%20.10e%20.10e%20.10e\n", (vx[i-1][j][k] + vx[i][j][k]) / 2.0, (vy[i][j-1][k] + vy[i][j][k]) / 2.0, (vz[i][j][k-1] + vz[i][j][k]) / 2.0);
}
}
}
}
}
fclose(fp);
}
//■進捗の出力■
if(t >= (int)((double)tcount * txstep))
{
printf("%s%7.2f%s\n", "Completed ", (double)t / (double)tx * 100.0, "%");
tcount++;
}
}
//■■■メモリの破棄■■■
free(p);
free(px);
free(py);
free(pz);
free(vx);
free(vy);
free(vz);
free(q);
return 0;
}
Raw
fdtd_3d_v2.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
print vx
# 音響領域
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])
print vx
# 左側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])
print vx
# -- (上) 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(j, 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(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
# 結果の出力
fcount += 1
filename = "%s_%05d.vtk" % (fn, fcount)
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%10d%s" % ("POINT_DATA", (ix - 2 * pl) * (jx - 2 * pl) * (k - 2 * pl), 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment