?? cpmlfdtd3d.c
字號:
/*******************************************************************
// 3-D FDTD code with CPML absorbing boundary conditions - Dipole
//******************************************************************
//
// This C code implements the finite-difference time-domain
// solution of Maxwell's curl equations over a three-dimensional
// Cartesian space lattice. The grid is terminated by CPML ABCs.
// The thickness of the PML in each Cartesian direction can be
// varied independently.
//
//******************************************************************/
//Header files (Libraries to be included)
#include<stdio.h>
#include<math.h>
#include<string.h>
#include<stdlib.h>
// Fundamental Constants (MKS units)
double pi = 3.14159265358979;
double C = 2.99792458E8;
double muO;
double epsO;
// Specify Material Relative Permittivity and Conductivity
double epsR = 1.0;//free space
// Specify Number of Time Steps and Grid Size Parameters
int nMax = 500; // total number of time steps
// grid size corresponding to the number of Ez field components
int Imax = 51;
int Jmax = 126;
int Kmax = 26;
// Specify Grid Cell Size in Each Direction and Calculate the
// Resulting Courant-Stable Time Step
double dx = 1.0E-3;
double dy = 1.0E-3;
double dz = 1.0E-3; // cell size in each direction
// time step increment
double dt;
// Specify the Impulsive Source (Differentiated Gaussian) parameters
double tw = 53.0E-12;//pulse width
double tO;//delay
double source;//Differentiated Gaussian source
double amp = 1000;// Amplitude
//Specify the Time Step at which the data has to be saved for Visualization
int save_modulus = 10;
// Specify the dipole Boundaries(A cuboidal rode- NOT as a cylinder)
int istart, iend, jstart;
int jend, kstart, kend;
//Output recording point
int ksource = 12;
// Specify the CPML Thickness in Each Direction (Value of Zero
// Corresponds to No PML, and the Grid is Terminated with a PEC)
// PML thickness in each direction
int nxPML_1, nxPML_2, nyPML_1;
int nyPML_2, nzPML_1, nzPML_2;
// Specify the CPML Order and Other Parameters:
int m = 3, ma = 1;
double sig_x_max;
double sig_y_max;
double sig_z_max;
double alpha_x_max;
double alpha_y_max;
double alpha_z_max;
double kappa_x_max;
double kappa_y_max;
double kappa_z_max;
//Loop indices
int i,j,ii,jj,k,kk,n;
// H & E Field components
float ***Hx;
float ***Hy;
float ***Hz;
float ***Ex;
float ***Ey;
float ***Ez;
short ***ID1;//medium definition array for Ex
short ***ID2;//medium definition array for Ey
short ***ID3;//medium definition array for Ez
// CPML components (Taflove 3rd Edition, Chapter 7)
float ***psi_Ezx_1;
float ***psi_Ezx_2;
float ***psi_Hyx_1;
float ***psi_Hyx_2;
float ***psi_Ezy_1;
float ***psi_Ezy_2;
float ***psi_Hxy_1;
float ***psi_Hxy_2;
float ***psi_Hxz_1;
float ***psi_Hxz_2;
float ***psi_Hyz_1;
float ***psi_Hyz_2;
float ***psi_Exz_1;
float ***psi_Exz_2;
float ***psi_Eyz_1;
float ***psi_Eyz_2;
float ***psi_Hzx_1;
float ***psi_Eyx_1;
float ***psi_Hzx_2;
float ***psi_Eyx_2;
float ***psi_Hzy_1;
float ***psi_Exy_1;
float ***psi_Hzy_2;
float ***psi_Exy_2;
float *be_x_1, *ce_x_1, *alphae_x_PML_1, *sige_x_PML_1, *kappae_x_PML_1;
float *bh_x_1, *ch_x_1, *alphah_x_PML_1, *sigh_x_PML_1, *kappah_x_PML_1;
float *be_x_2, *ce_x_2, *alphae_x_PML_2, *sige_x_PML_2, *kappae_x_PML_2;
float *bh_x_2, *ch_x_2, *alphah_x_PML_2, *sigh_x_PML_2, *kappah_x_PML_2;
float *be_y_1, *ce_y_1, *alphae_y_PML_1, *sige_y_PML_1, *kappae_y_PML_1;
float *bh_y_1, *ch_y_1, *alphah_y_PML_1, *sigh_y_PML_1, *kappah_y_PML_1;
float *be_y_2, *ce_y_2, *alphae_y_PML_2, *sige_y_PML_2, *kappae_y_PML_2;
float *bh_y_2, *ch_y_2, *alphah_y_PML_2, *sigh_y_PML_2, *kappah_y_PML_2;
float *be_z_1, *ce_z_1, *alphae_z_PML_1, *sige_z_PML_1, *kappae_z_PML_1;
float *bh_z_1, *ch_z_1, *alphah_z_PML_1, *sigh_z_PML_1, *kappah_z_PML_1;
float *be_z_2, *ce_z_2, *alphae_z_PML_2, *sige_z_PML_2, *kappae_z_PML_2;
float *bh_z_2, *ch_z_2, *alphah_z_PML_2, *sigh_z_PML_2, *kappah_z_PML_2;
// denominators for the update equations
float *den_ex;
float *den_hx;
float *den_ey;
float *den_hy;
float *den_ez;
float *den_hz;
//Max number of materials allowed
int numMaterials = 50;
//permittivity, permeability and conductivity of diffrent materials
double *epsilon;
double *mu;
double *sigma;
//E field update coefficients
float *CA;
float *CB;
//H field update coefficients
float DA;
float DB;
//Function prototype definitions
void initialize();//Memeory initialization
void setUp();//Coefficients, parameters etc will get computed
void initializeCPML();//CPML coefficient computation
void compute();//E & H Field update equation
void buildObject();//Creates the object geometry
void yeeCube (int, int,int, short);//Sets material properties to a cell
void writeField(int);//Writes output
void buildSphere();//Builds a spherical object
void buildDipole();//Builds a dipole
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void main() {
initialize();
setUp();
buildObject();
initializeCPML();
compute();
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void initialize() {
muO = 4.0 * pi * 1.0E-7;
epsO = 1.0 / (C * C * muO);
//PML Layers (10 layers)
nxPML_1 = 11;
nxPML_2 = 11;
nyPML_1 = 11;
nyPML_2 = 11;
nzPML_1 = 11;
nzPML_2 = 11;
//Dynamic memory allocation
epsilon = (double *)malloc((numMaterials) * sizeof(double));
for(i = 0; i < numMaterials; i++) {
epsilon[i] = epsO;
}
mu = (double *)malloc((numMaterials) * sizeof(double));
for(i = 0; i < numMaterials; i++) {
mu[i] = muO;
}
sigma = (double *)malloc((numMaterials) * sizeof(double));
for(i = 0; i < numMaterials; i++) {
sigma[i] = 0.0;
}
CA = (float *)malloc((numMaterials) * sizeof(float));
for(i = 0; i < numMaterials; i++) {
CA[i] = 0.0;
}
CB = (float *)malloc((numMaterials) * sizeof(float));
for(i = 0; i < numMaterials; i++) {
CB[i] = 0.0;
}
Ez = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
Ez[i] = (float **)malloc(Jmax * sizeof(float *));
for(j = 0; j < Jmax; j++) {
Ez[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
Ez[i][j][k] = 0.0;
}
}
}
Ey = (float ***)malloc((Imax) * sizeof(float **));
for(i = 0; i < Imax; i++) {
Ey[i] = (float **)malloc((Jmax-1) * sizeof(float *));
for(j = 0; j < Jmax -1; j++) {
Ey[i][j] = (float *)malloc((Kmax-1) * sizeof(float));
for(k = 0; k < Kmax-1; k++) {
Ey[i][j][k] = 0.0;
}
}
}
Ex = (float ***)malloc((Imax-1) * sizeof(float **));
for(i = 0; i < Imax-1; i++) {
Ex[i] = (float **)malloc(Jmax * sizeof(float *));
for(j = 0; j < Jmax; j++) {
Ex[i][j] = (float *)malloc((Kmax-1) * sizeof(float));
for(k = 0; k < Kmax-1; k++) {
Ex[i][j][k] = 0.0;
}
}
}
Hx = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
Hx[i] = (float **)malloc((Jmax-1) * sizeof(float *));
for(j = 0; j < Jmax-1; j++) {
Hx[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
Hx[i][j][k] = 0.0;
}
}
}
Hy = (float ***)malloc((Imax-1) * sizeof(float **));
for(i = 0; i < Imax-1; i++) {
Hy[i] = (float **)malloc(Jmax * sizeof(float *));
for(j = 0; j < Jmax; j++) {
Hy[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
Hy[i][j][k] = 0.0;
}
}
}
Hz = (float ***)malloc((Imax-1) * sizeof(float **));
for(i = 0; i < Imax-1; i++) {
Hz[i] = (float **)malloc((Jmax-1) * sizeof(float *));
for(j = 0; j < Jmax-1; j++) {
Hz[i][j] = (float *)malloc((Kmax-1) * sizeof(float));
for(k = 0; k < Kmax-1; k++) {
Hz[i][j][k] = 0.0;
}
}
}
ID1 = (short ***)malloc(Imax * sizeof(short **));
for(i = 0; i < Imax; i++) {
ID1[i] = (short **)malloc(Jmax * sizeof(short *));
for(j = 0; j < Jmax; j++) {
ID1[i][j] = (short *)malloc(Kmax * sizeof(short));
for(k = 0; k < Kmax; k++) {
ID1[i][j][k] = 0;
}
}
}
ID2 = (short ***)malloc(Imax * sizeof(short **));
for(i = 0; i < Imax; i++) {
ID2[i] = (short **)malloc(Jmax * sizeof(short *));
for(j = 0; j < Jmax; j++) {
ID2[i][j] = (short *)malloc(Kmax * sizeof(short));
for(k = 0; k < Kmax; k++) {
ID2[i][j][k] = 0;
}
}
}
ID3 = (short ***)malloc(Imax * sizeof(short **));
for(i = 0; i < Imax; i++) {
ID3[i] = (short **)malloc(Jmax * sizeof(short *));
for(j = 0; j < Jmax; j++) {
ID3[i][j] = (short *)malloc(Kmax * sizeof(short));
for(k = 0; k < Kmax; k++) {
ID3[i][j][k] = 0;
}
}
}
psi_Ezx_1 = (float ***)malloc(nxPML_1 * sizeof(float **));
for(i = 0; i < nxPML_1; i++) {
psi_Ezx_1[i] = (float **)malloc(Jmax * sizeof(float *));
for(j = 0; j < Jmax; j++) {
psi_Ezx_1[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Ezx_1[i][j][k] = 0.0;
}
}
}
psi_Ezx_2 = (float ***)malloc(nxPML_2 * sizeof(float **));
for(i = 0; i < nxPML_2; i++) {
psi_Ezx_2[i] = (float **)malloc(Jmax * sizeof(float *));
for(j = 0; j < Jmax; j++) {
psi_Ezx_2[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Ezx_2[i][j][k] = 0.0;
}
}
}
psi_Hyx_1 = (float ***)malloc((nxPML_1-1) * sizeof(float **));
for(i = 0; i < nxPML_1-1; i++) {
psi_Hyx_1[i] = (float **)malloc(Jmax * sizeof(float *));
for(j = 0; j < Jmax; j++) {
psi_Hyx_1[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Hyx_1[i][j][k] = 0.0;
}
}
}
psi_Hyx_2 = (float ***)malloc((nxPML_2-1) * sizeof(float **));
for(i = 0; i < nxPML_1-1; i++) {
psi_Hyx_2[i] = (float **)malloc(Jmax * sizeof(float *));
for(j = 0; j < Jmax; j++) {
psi_Hyx_2[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Hyx_2[i][j][k] = 0.0;
}
}
}
psi_Ezy_1 = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
psi_Ezy_1[i] = (float **)malloc(nyPML_1 * sizeof(float *));
for(j = 0; j < nyPML_1; j++) {
psi_Ezy_1[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Ezy_1[i][j][k] = 0.0;
}
}
}
psi_Ezy_2 = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
psi_Ezy_2[i] = (float **)malloc(nyPML_2 * sizeof(float *));
for(j = 0; j < nyPML_2; j++) {
psi_Ezy_2[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Ezy_2[i][j][k] = 0.0;
}
}
}
psi_Hxy_1 = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
psi_Hxy_1[i] = (float **)malloc((nyPML_1-1) * sizeof(float *));
for(j = 0; j < nyPML_1-1; j++) {
psi_Hxy_1[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Hxy_1[i][j][k] = 0.0;
}
}
}
psi_Hxy_2 = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
psi_Hxy_2[i] = (float **)malloc((nyPML_2-1) * sizeof(float *));
for(j = 0; j < nyPML_2-1; j++) {
psi_Hxy_2[i][j] = (float *)malloc(Kmax * sizeof(float));
for(k = 0; k < Kmax; k++) {
psi_Hxy_2[i][j][k] = 0.0;
}
}
}
psi_Hxz_1 = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
psi_Hxz_1[i] = (float **)malloc((Jmax-1) * sizeof(float *));
for(j = 0; j < Jmax; j++) {
psi_Hxz_1[i][j] = (float *)malloc((nzPML_1-1) * sizeof(float));
for(k = 0; k < nzPML_1-1; k++) {
psi_Hxz_1[i][j][k] = 0.0;
}
}
}
psi_Hxz_2 = (float ***)malloc(Imax * sizeof(float **));
for(i = 0; i < Imax; i++) {
psi_Hxz_2[i] = (float **)malloc((Jmax-1) * sizeof(float *));
for(j = 0; j < Jmax; j++) {
psi_Hxz_2[i][j] = (float *)malloc((nzPML_2-1) * sizeof(float));
for(k = 0; k < nzPML_2-1; k++) {
psi_Hxz_2[i][j][k] = 0.0;
}
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -