sixfu
Newbie level 1
2d pml
I am learning FDTD through Sullivan's book---Electromagnetic Simulation Using the FDTD Method, and I can not get the right result from the code listed in the book for fd2d_3.2.c. Also I downloaded this code from his website, it doesn't work too. Can any body help me to figure it out!
Thanks a lot!
Below is the source code from the book:
*******************************************************
/* Fd2d_3.2.c. 2D TM program with the PML */
# include <math.h>
# include <stdlib.h>
# include <stdio.h>
#define IE 140
#define JE 140
main ()
{
float ga[IE][JE],dz[IE][JE],ez[IE][JE],hx[IE][JE],hy[IE][JE];
int l,n,i,j,ic,jc,nsteps,npml;
float ddx,dt,T,epsz,pi,epsilon,sigma,eaf;
float xn,xxn,xnum,xd,curl_e;
float t0,spread,pulse;
float gi2[IE],gi3[IE];
float gj2[JE],gj3[IE];
float fi1[IE],fi2[IE],fi3[JE];
float fj1[JE],fj2[JE],fj3[JE];
float ihx[IE][JE],ihy[IE][JE];
FILE *fp, *fopen();
ic = IE/2-20;
jc = JE/2-20;
ddx = .01; /* Cell size */
dt =ddx/6e8; /* Time steps */
epsz = 8.8e-12;
pi=3.14159;
/* Initialize the arrays */
for ( j=0; j < JE; j++ ) {
printf( "%2d ",j);
for ( i=0; i < IE; i++ ) {
dz[j]= 0.0 ;
hx[j]= 0.0 ;
hy[j]= 0.0 ;
ihx[j]= 0.0 ;
ihy[j]= 0.0 ;
ga[j]= 1.0 ;
printf( "%5.2f ",ga[j]);
}
printf( " \n");
}
/* Calculate the PML parameters */
for (i=0;i< IE; i++) {
gi2 = 1.0;
gi3 = 1.0;
fi1 = 0.0;
fi2 = 1.0;
fi3 = 1.0;
}
for (j=0;j< IE; j++) {
gj2[j] = 1.0;
gj3[j] = 1.0;
fj1[j] = 0.0;
fj2[j] = 1.0;
fj3[j] = 1.0;
}
printf( "Number of PML cells --> ");
scanf("%d", &npml);
for (i=0;i<= npml; i++) {
xnum = npml - i;
xd = npml;
xxn = xnum/xd;
xn = 0.25*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
gi2 = 1.0/(1.0+xn);
gi2[IE-1-i] = 1.0/(1.0+xn);
gi3 = (1.0 - xn)/(1.0 + xn);
gi3[IE-i-1] = (1.0 - xn)/(1.0 + xn);
xxn = (xnum-.5)/xd;
xn = 0.25*pow(xxn,3.0);
fi1 = xn;
fi1[IE-2-i] = xn;
fi2 = 1.0/(1.0+xn);
fi2[IE-2-i] = 1.0/(1.0+xn);
fi3 = (1.0 - xn)/(1.0 + xn);
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
}
for (j=0;j<= npml; j++) {
xnum = npml - j;
xd = npml;
xxn = xnum/xd;
xn = 0.25*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
gj2[j] = 1.0/(1.0+xn);
gj2[JE-1-j] = 1.0/(1.0+xn);
gj3[j] = (1.0 - xn)/(1.0 + xn);
gj3[JE-j-1] = (1.0 - xn)/(1.0 + xn);
xxn = (xnum-.5)/xd;
xn = 0.25*pow(xxn,3.0);
fj1[j] = xn;
fj1[JE-2-j] = xn;
fj2[j] = 1.0/(1.0+xn);
fj2[JE-2-j] = 1.0/(1.0+xn);
fj3[j] = (1.0 - xn)/(1.0 + xn);
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
}
printf("gi + fi \n");
for (i=0; i< IE; i++) {
printf( "%2d %5.2f %5.2f \n",
i,gi2,gi3),
printf( " %5.2f %5.2f %5.2f \n ",
fi1,fi2,fi3);
}
printf("gj + fj \n");
for (j=0; j< JE; j++) {
printf( "%2d %5.2f %5.2f \n",
j,gj2[j],gj3[j]),
printf( " %5.2f %5.2f %5.2f \n ",
fj1[j],fj2[j],fj3[j]);
}
t0 = 40.0;
spread = 12.0;
T = 0;
nsteps = 1;
while ( nsteps > 0 ) {
printf( "nsteps --> ");
scanf("%d", &nsteps);
printf("%d \n", nsteps);
for ( n=1; n <=nsteps ; n++) {
T = T + 1;
/* ---- Start of the Main FDTD loop ---- */
/* Calculate the Dz field */
for ( j=1; j < IE; j++ ) {
for ( i=1; i < IE; i++ ) {
dz[j] = gi3*gj3[j]*dz[j]
+ gi2*gj2[j]*.5*( hy[j] - hy[i-1][j]
- hx[j] + hx[j-1]) ;
}
}
/* Sinusoidal Source */
/* pulse = sin(2*pi*1500*1e6*dt*T);; */
pulse = exp(-.5*pow( (T-t0)/spread,2.));
dz[ic][jc] = pulse;
/* Calculate the Ez field */
/* Leave the Ez edges to 0, as part of the PML */
for ( j=1; j < JE-1; j++ ) {
for ( i=1; i < IE-1; i++ ) {
ez[j] = ga[j]*dz[j] ;
}
}
printf("%3f %6.2f \n ",T,ez[ic][jc]);
/* Calculate the Hx field */
for ( j=0; j < JE-1; j++ ) {
for ( i=0; i < IE; i++ ) {
curl_e = ez[j] - ez[j+1] ;
ihx[j] = ihx[j] + fi1*curl_e ;
hx[j] = fj3[j]*hx[j]
+ fj2[j]*.5*(curl_e + ihx[j]);
}
}
/* Calculate the Hy field */
for ( j=0; j <= JE-1; j++ ) {
for ( i=0; i < IE-1; i++ ) {
curl_e = ez[i+1][j] - ez[j];
ihy[j] = ihy[j] + fj1[j]*curl_e ;
hy[j] = fi3*hy[j]
+ fi2*.5*(curl_e + ihy[j]);
}
}
}
/* ---- End of the main FDTD loop ---- */
for ( j=1; j < JE; j++ ) {
printf( "%2d ",j);
for ( i=1; i <= IE; i++ ) {
printf( "%4.1f",ez[j]);
}
printf( " \n");
}
/* Write the E field out to a file "Ez" */
fp = fopen( "Ez","w");
for ( j=0; j < JE; j++ ) {
for ( i=0; i < IE; i++ ) {
fprintf( fp,"%6.3f ",ez[j]);
}
fprintf( fp," \n");
}
fclose(fp);
printf("T = %6.0f \n ",T);
}
}
I am learning FDTD through Sullivan's book---Electromagnetic Simulation Using the FDTD Method, and I can not get the right result from the code listed in the book for fd2d_3.2.c. Also I downloaded this code from his website, it doesn't work too. Can any body help me to figure it out!
Thanks a lot!
Below is the source code from the book:
*******************************************************
/* Fd2d_3.2.c. 2D TM program with the PML */
# include <math.h>
# include <stdlib.h>
# include <stdio.h>
#define IE 140
#define JE 140
main ()
{
float ga[IE][JE],dz[IE][JE],ez[IE][JE],hx[IE][JE],hy[IE][JE];
int l,n,i,j,ic,jc,nsteps,npml;
float ddx,dt,T,epsz,pi,epsilon,sigma,eaf;
float xn,xxn,xnum,xd,curl_e;
float t0,spread,pulse;
float gi2[IE],gi3[IE];
float gj2[JE],gj3[IE];
float fi1[IE],fi2[IE],fi3[JE];
float fj1[JE],fj2[JE],fj3[JE];
float ihx[IE][JE],ihy[IE][JE];
FILE *fp, *fopen();
ic = IE/2-20;
jc = JE/2-20;
ddx = .01; /* Cell size */
dt =ddx/6e8; /* Time steps */
epsz = 8.8e-12;
pi=3.14159;
/* Initialize the arrays */
for ( j=0; j < JE; j++ ) {
printf( "%2d ",j);
for ( i=0; i < IE; i++ ) {
dz[j]= 0.0 ;
hx[j]= 0.0 ;
hy[j]= 0.0 ;
ihx[j]= 0.0 ;
ihy[j]= 0.0 ;
ga[j]= 1.0 ;
printf( "%5.2f ",ga[j]);
}
printf( " \n");
}
/* Calculate the PML parameters */
for (i=0;i< IE; i++) {
gi2 = 1.0;
gi3 = 1.0;
fi1 = 0.0;
fi2 = 1.0;
fi3 = 1.0;
}
for (j=0;j< IE; j++) {
gj2[j] = 1.0;
gj3[j] = 1.0;
fj1[j] = 0.0;
fj2[j] = 1.0;
fj3[j] = 1.0;
}
printf( "Number of PML cells --> ");
scanf("%d", &npml);
for (i=0;i<= npml; i++) {
xnum = npml - i;
xd = npml;
xxn = xnum/xd;
xn = 0.25*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
gi2 = 1.0/(1.0+xn);
gi2[IE-1-i] = 1.0/(1.0+xn);
gi3 = (1.0 - xn)/(1.0 + xn);
gi3[IE-i-1] = (1.0 - xn)/(1.0 + xn);
xxn = (xnum-.5)/xd;
xn = 0.25*pow(xxn,3.0);
fi1 = xn;
fi1[IE-2-i] = xn;
fi2 = 1.0/(1.0+xn);
fi2[IE-2-i] = 1.0/(1.0+xn);
fi3 = (1.0 - xn)/(1.0 + xn);
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
}
for (j=0;j<= npml; j++) {
xnum = npml - j;
xd = npml;
xxn = xnum/xd;
xn = 0.25*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
gj2[j] = 1.0/(1.0+xn);
gj2[JE-1-j] = 1.0/(1.0+xn);
gj3[j] = (1.0 - xn)/(1.0 + xn);
gj3[JE-j-1] = (1.0 - xn)/(1.0 + xn);
xxn = (xnum-.5)/xd;
xn = 0.25*pow(xxn,3.0);
fj1[j] = xn;
fj1[JE-2-j] = xn;
fj2[j] = 1.0/(1.0+xn);
fj2[JE-2-j] = 1.0/(1.0+xn);
fj3[j] = (1.0 - xn)/(1.0 + xn);
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
}
printf("gi + fi \n");
for (i=0; i< IE; i++) {
printf( "%2d %5.2f %5.2f \n",
i,gi2,gi3),
printf( " %5.2f %5.2f %5.2f \n ",
fi1,fi2,fi3);
}
printf("gj + fj \n");
for (j=0; j< JE; j++) {
printf( "%2d %5.2f %5.2f \n",
j,gj2[j],gj3[j]),
printf( " %5.2f %5.2f %5.2f \n ",
fj1[j],fj2[j],fj3[j]);
}
t0 = 40.0;
spread = 12.0;
T = 0;
nsteps = 1;
while ( nsteps > 0 ) {
printf( "nsteps --> ");
scanf("%d", &nsteps);
printf("%d \n", nsteps);
for ( n=1; n <=nsteps ; n++) {
T = T + 1;
/* ---- Start of the Main FDTD loop ---- */
/* Calculate the Dz field */
for ( j=1; j < IE; j++ ) {
for ( i=1; i < IE; i++ ) {
dz[j] = gi3*gj3[j]*dz[j]
+ gi2*gj2[j]*.5*( hy[j] - hy[i-1][j]
- hx[j] + hx[j-1]) ;
}
}
/* Sinusoidal Source */
/* pulse = sin(2*pi*1500*1e6*dt*T);; */
pulse = exp(-.5*pow( (T-t0)/spread,2.));
dz[ic][jc] = pulse;
/* Calculate the Ez field */
/* Leave the Ez edges to 0, as part of the PML */
for ( j=1; j < JE-1; j++ ) {
for ( i=1; i < IE-1; i++ ) {
ez[j] = ga[j]*dz[j] ;
}
}
printf("%3f %6.2f \n ",T,ez[ic][jc]);
/* Calculate the Hx field */
for ( j=0; j < JE-1; j++ ) {
for ( i=0; i < IE; i++ ) {
curl_e = ez[j] - ez[j+1] ;
ihx[j] = ihx[j] + fi1*curl_e ;
hx[j] = fj3[j]*hx[j]
+ fj2[j]*.5*(curl_e + ihx[j]);
}
}
/* Calculate the Hy field */
for ( j=0; j <= JE-1; j++ ) {
for ( i=0; i < IE-1; i++ ) {
curl_e = ez[i+1][j] - ez[j];
ihy[j] = ihy[j] + fj1[j]*curl_e ;
hy[j] = fi3*hy[j]
+ fi2*.5*(curl_e + ihy[j]);
}
}
}
/* ---- End of the main FDTD loop ---- */
for ( j=1; j < JE; j++ ) {
printf( "%2d ",j);
for ( i=1; i <= IE; i++ ) {
printf( "%4.1f",ez[j]);
}
printf( " \n");
}
/* Write the E field out to a file "Ez" */
fp = fopen( "Ez","w");
for ( j=0; j < JE; j++ ) {
for ( i=0; i < IE; i++ ) {
fprintf( fp,"%6.3f ",ez[j]);
}
fprintf( fp," \n");
}
fclose(fp);
printf("T = %6.0f \n ",T);
}
}