Continue to Site

Welcome to EDAboard.com

Welcome to our site! EDAboard.com is an international Electronics Discussion Forum focused on EDA software, circuits, schematics, books, theory, papers, asic, pld, 8051, DSP, Network, RF, Analog Design, PCB, Service Manuals... and a whole lot more! To participate you need to register. Registration is free. Click here to register now.

Register Log in

Problems from Sullivan's FDTD 3D codes

Status
Not open for further replies.
lqkhai

lqkhai

Junior Member level 3
Joined
Mar 7, 2006
Messages
31
Helped
3
Reputation
6
Reaction score
0
Trophy points
1,286
Activity points
1,465
Dear,
I met some problems with FDTD 3D codes.
This is my codes.
Code: 
/*File: FD3D.h*/

#include "iostream.h"
#include "fstream.h"
#include "math.h"
#include "stdlib.h"
#include "malloc.h"

#define IE 40
#define JE 40
#define KE 40
#define N 40


class ThreeDFDTD
{
	public: 
		/*float gax[IE][JE][KE],gay[IE][JE][KE],gaz[IE][JE][KE];
		float dx[IE][JE][KE],dy[IE][JE][KE],dz[IE][JE][KE];
		float ex[IE][JE][KE],ey[IE][JE][KE],ez[IE][JE][KE];
		float hx[IE][JE][KE],hy[IE][JE][KE],hz[IE][JE][KE];*/
		float ***gax,***gay,***gaz;
		float ***dx,***dy,***dz;
		float ***ex,***ey,***ez;
		float ***hx,***hy,***hz;
		float ddx,dt,T,epsz,epsilon,sigma,eaf;
		float xn,xxn,xnum,xd,curl_e;
		float t0,spread,pulse,NSTEPs;
		int ic,jc,kc;
		void createMang(float ***&mang);
		void freeMang(float ***&mang);
		void initialize(void);
		void fdtd(void);
};

Code: 
/*File: Main.cpp*/

#include "FD3D.h"

int main()
{
	int mode;
	ThreeDFDTD c;

	while(1)
	{
		cout << endl
			 << "****3 Dimensional FDTD simulation*****" 
			 << "\n****1> 3D FDTD. Dipole in free space"   
			 << "\n****2> Exit" 
		     << endl;
		cin  >> mode;
		switch(mode)
		{
		case 1:
			{	
			c.fdtd();
			break;
			}
		default:exit(0);
		}
	}
	return 0;
}

Code: 
/*File FD3D.cpp*/

#include "FD3D.h"

void ThreeDFDTD::createMang(float ***&mang)
{
      /* mang = new double**[N];
for (int i=0; i<N; i++)
    mang[i] = new double*[N];

for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
    mang[i][j] = new double[N];*/
// gan gia tri
/*for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
for (int k=0; k<N; k++)
    mang[i][j][k] = i+j+k;*/
    /*mang=(float***)calloc(IE*JE*KE,sizeof(float**));
    for(int i=0;i<N;i++)
    {
            mang[i]=(float**)calloc(IE*JE,sizeof(float*));
            for(int j=0;j<N;j++)
            {
                 mang[i][j]=(float*)calloc(IE,sizeof(float));
            }
    }*/
    mang=(float***)malloc(IE*JE*KE*sizeof(float**));
    for(int i=0;i<N;i++)
    {
            mang[i]=(float**)malloc(IE*JE*sizeof(float*));
            for(int j=0;j<N;j++)
            {
                 mang[i][j]=(float*)malloc(IE*sizeof(float));
            }
    }   
}

void ThreeDFDTD::freeMang(float ***&mang)
{
// giai phong mang
/*for (int i=0; i<N; i++)
for (int j=0; j<N; j++)
    delete mang[i][j];
    
for (int i=0; i<N; i++)
    delete mang[i];

     delete mang;*/
     for(int i=0;i<N;i++)
     {
          for(int j=0;j<N;j++)
              free(mang[i][j]);
          free(mang[i]);
     }
     free(mang);
     
}

void ThreeDFDTD::initialize(void)
{
	int i,j,k;
	ic=IE/2;
	jc=JE/2;
	kc=KE/2;
	ddx=.01;
	dt=ddx/(6e8);
	epsz=8.8e-12;
	/*createMang(ex);
	createMang(ey);
	createMang(ez);
	createMang(hx);
	createMang(hy);
	createMang(hz);
	createMang(dx);
	createMang(dy);
	createMang(dz);
	createMang(gax);
	createMang(gay);
	createMang(gaz);
	
	
	for(k=0;k<KE;k++)
		for(j=0;j<JE;j++)
			for(i=0;i<IE;i++)
			{
				ex[i][j][k]=0.0;
				ey[i][j][k]=0.0;
				ez[i][j][k]=0.0;
				dx[i][j][k]=0.0;
				dy[i][j][k]=0.0;
				dz[i][j][k]=0.0;
				hx[i][j][k]=0.0;
				hy[i][j][k]=0.0;
				hz[i][j][k]=0.0;
				gax[i][j][k]=1.0;
				gay[i][j][k]=1.0;
				gaz[i][j][k]=1.0;
			}

	for(k=11;k<30;k++)
	{
		gaz[ic][jc][k]=0.;
	}
	gaz[ic][jc][kc]=0.;
	t0=20.0;
	spread=6.0;
	T=0;
	NSTEPs=1;*/
}
void ThreeDFDTD::fdtd(void)
{
	int n;
	int i,j,k;
	initialize();
	
	createMang(ex);
	createMang(ey);
	createMang(ez);
	createMang(hx);
	createMang(hy);
	createMang(hz);
	createMang(dx);
	createMang(dy);
	createMang(dz);
	createMang(gax);
	createMang(gay);
	createMang(gaz);
	
	
	for(k=0;k<KE;k++)
		for(j=0;j<JE;j++)
			for(i=0;i<IE;i++)
			{
				ex[i][j][k]=0.0;
				ey[i][j][k]=0.0;
				ez[i][j][k]=0.0;
				dx[i][j][k]=0.0;
				dy[i][j][k]=0.0;
				dz[i][j][k]=0.0;
				hx[i][j][k]=0.0;
				hy[i][j][k]=0.0;
				hz[i][j][k]=0.0;
				gax[i][j][k]=1.0;
				gay[i][j][k]=1.0;
				gaz[i][j][k]=1.0;
			}

	for(k=11;k<30;k++)
	{
		gaz[ic][jc][k]=0.;
	}
	//gaz[ic][jc][kc]=0.;
	t0=20.0;
	spread=6.0;
	T=0;
	NSTEPs=1;
	
	cout << endl
			 << "Please input NSTEPs= ";
		cin  >> NSTEPs;
		
		for(n=1;n<=NSTEPs;n++)
		{
			T=T+1;
			/****Start of the main FDTD loop****/
			/****Calculate the Dx field****/
			for(k=1;k<KE;k++)
				for(j=1;j<JE;j++)
					for(i=1;i<IE;i++)
					{
					dx[i][j][k]=dx[i][j][k]+.5*(hz[i][j][k]-hz[i][j-1][k]-hy[i][j][k]+hy[i][j][k-1]);
					}
			/****Calculate the Dy field****/
			for(k=1;k<KE;k++)
				for(j=1;j<JE;j++)
					for(i=1;i<IE;i++)
					{
					dy[i][j][k]=dy[i][j][k]+.5*(hx[i][j][k]-hx[i][j][k-1]-hz[i][j][k]+hz[i-1][j][k]);
					}

			/****Calculate the Dz field****/
			for(k=1;k<KE;k++)
				for(j=1;j<JE;j++)
					for(i=1;i<IE;i++)
					{
					dz[i][j][k]=dz[i][j][k]+.5*(hy[i][j][k]-hy[i-1][j][k]-hx[i][j][k]+hx[i][j-1][k]);
					}
			/****Source****/
			pulse=exp(-.5*pow((t0-T)/spread,2));
			dz[ic][jc][kc]=pulse;
			/****Calculate the E field from D****/
			for(k=1;k<KE-1;k++)
				for(j=1;j<JE-1;j++)
					for(i=1;i<IE-1;i++)
					{
						ex[i][j][k]=gax[i][j][k]*dx[i][j][k];
						ey[i][j][k]=gay[i][j][k]*dy[i][j][k];
						ez[i][j][k]=gaz[i][j][k]*dz[i][j][k];
					}
			/****Calculate the Hx field****/
			for(k=1;k<KE-1;k++)
				for(j=1;j<JE-1;j++)
					for(i=1;i<IE;i++)
					{
					hx[i][j][k]=hx[i][j][k]+.5*(ey[i][j][k+1]-ey[i][j][k]-ez[i][j+1][k]-ez[i][j][k]);
					}
			/****Calculate the Hy field****/
			for(k=1;k<KE-1;k++)
				for(j=1;j<JE;j++)
					for(i=1;i<IE-1;i++)
					{
					hy[i][j][k]=hy[i][j][k]+.5*(ez[i+1][j][k]-ez[i][j][k]-ex[i][j][k+1]+ex[i][j][k]);
					}
			/****Calculate the Hx field****/
			for(k=1;k<KE;k++)
				for(j=1;j<JE-1;j++)
					for(i=1;i<IE-1;i++)
					{
					hz[i][j][k]=hz[i][j][k]+.5*(ex[i][j+1][k]-ex[i][j][k]-ey[i+1][j][k]+ey[i][j][k]);
					}
			/****End of the main FDTD loop****/
			/****Write the Ez field into Ez.txt****/
					ofstream Ez("Ez.txt");
			for(j=0;j<JE;j++)
			{
				for(i=0;i<IE;i++)
				{
					 Ez <<   ez[i][j][kc] <<" ";
				}
				Ez << "\n";
			}
			Ez.close();
		
			for(k=1;k<KE;k++)
				for(i=1;i<IE-1;i++)
				{
					cout << ez[i][jc][k] << " ";
				}
				cout << endl;
		}
		
    printf("\n No error there");
    
	freeMang(ex);
	freeMang(ey);
	freeMang(ez);
	freeMang(hx);
	freeMang(hy);
	freeMang(hz);
	freeMang(dx);
	freeMang(dy);
	freeMang(dz);
	freeMang(gax);
	freeMang(gay);
	freeMang(gaz);
	
}

In my codes, If I use static 3D arrays as Sullivan's codes, The error will occur as the following: Array size too large.
In order to deal with memory problems, I have to use the Dynamic 3D arrays as my codes. But if the array size is 20x20x20, it's ok. In the case of array size of 40x40x40, the results are all 0.0.

Please help me.

cheers,

lqkhai
 

What is the size of your ram increase it and try again and be careful about not using as for example ex[IE][JE][KE] this will be wrong sometimes because of this you can receive error by the way if you have other codes written please can give me thanks lets get in touch because i also try to follow the same book for my thesis

Added after 5 minutes:

I think you don't have such mistake that i mentioned but if i notice i will inform you see you soon now i am studying pml in 2D and then need to apply it for 3D and do my simulations of conducting wires by the way your programming is very good i think and you try to convert the normal programs in hte book to the function form , don't you have ever tried the program without functions , and then step by step write the as functions then you can cathc the mistake if you have , see you soon
 

Thanks for your attention!

don't you have ever tried the program without functions , and then step by step write the as functions then you can cathc the mistake if you have , see you soon
Click to expand...

I don't think whether my program with function form met problems. Anyway, if I code without function form the error also occurs as below.
I also try to convert 3D form into 1D as the following
Code: 
//Allocate dynamic memory
#define ALLOC_3D(name,nx,ny,nz,type){ \
        name=(type*)calloc(nx*ny*nz,sizeof(type));\
        if(!name){                        \
                 perror("ALLOC_3D");      \
                 exit(-1);                \
        };                                \
}                                         

#define Ez(M,N,P) ez[(M*IE+N)*JE+P]

int main()
{
   
    double *ez;
    ALLOC_3D(ez,IE,JE,KE,double);
  ......
}

But my final figure is different from Sullivan's figure 4.3(page 82).

cheers,

lqkhai
 

Status
Not open for further replies.

Similar threads

Part and Inventory Search

Welcome to EDABoard.com

Sponsor

Back
Top