fdtd_3d_xyzpml_threads_decomp.cpp

利用c++语言写的三维FDTD

/////////////////////////////////////////////////////////////
//set as follows:
//MALLOCMULTIHEAP=true - to allow the use of multiple heaps of free memory, 
//                       rather than just one - threads
/////////////////////////////////////////////////////////////
#include <pthread.h>  //for threads
#include <unistd.h>   //for getwd()
#include <sys/stat.h> //for mkdir()
#include <errno.h>    //for error handling

#include <stdio.h>
#include <math.h>
#include <ctime>
#include <malloc.h>

#include "Matrix.h"
#include "Load_Save_File_Data.h"
#include "fdtd_3d_alloc.h"
#include "fdtd_3d.h"
#include "Thread_Calc.h"
#include "Init_Thread_Calc_Data.h"
#include "FDTD_1D_EH_PML_LOSS.h"
#include "Alloc_Save_Load_Field_Thread.h"
#include "fdtd_3d_free.h"
#include "generate_report.h"

#include <iostream>
using std::cout;
using std::endl;
using std::fixed;

extern int errno;//error handling

//#define eval_W; //comment out if the energy is not evaluated

#include "read_input_file.h"

int main(int argc, char* argv[])
{
	///////////////////////////////////////////////////////
	//initialization to read from file the input data
	///////////////////////////////////////////////////////
	char *file_name = NULL;
	file_name = (char *) calloc(512,sizeof(char));
	if (!file_name)
	{
		cout << "Memory allocation problem - *file_name" << endl;
		return 1;
	}
		
	if (argc > 1)
	{
		strcpy(file_name,argv[1]);
	}
	else
	{
		strcpy(file_name,"Input_Data.txt");
	}
	
	Input_Data *Inp_D = NULL;
	Inp_D = (Input_Data *) calloc(1,sizeof(Input_Data));

	Inp_D->Path_Load_Workspace_Data = NULL;
	Inp_D->path_name_CoordPointSource = NULL;
	Inp_D->path_name_Index = NULL;
	Inp_D->path_name_MatParam = NULL;
	Inp_D->path_name_ParamPointSource = NULL;
	Inp_D->path_name_Wigner_Cell = NULL;
    
	if (readFile_fdtd_data(Inp_D,file_name))
	{
		cout << "Error reading Input_Data.dat" << endl;
		return 1;
	}
	else
	{
		cout << "Input file: " << file_name << " loaded succesfull" << endl;
	}
	///////////////////////////////////////////////////////

	int rc; //for error handling
	
	//to measure the program performance
	//for Unix
	double time_start, time_end, el_time;
	struct timeval tv;
	struct timezone tz;

	gettimeofday(&tv, &tz);
	time_start = (double) tv.tv_sec + (double) tv.tv_usec / 1000000.0;

	errno = 0;   //set errno to  zero

	///////////////////////////////////////////////////////////////////////////////////////
	//some usefull constants
	///////////////////////////////////////////////////////////////////////////////////////
	double pi = 180.0*atan(1.0)/45.0;
	//permittivity of free space 
	double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    double mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]
	//the speed of the light in the vacuum
	double speed_c = 1.0/sqrt(eps_0*mu_0); // [m/s]
	///////////////////////////////////////////////////////////////////////////////////////

	///////////////////////////////////////////////////////////////////////////////////////
	//the parameters of the photonic crystal
	///////////////////////////////////////////////////////////////////////////////////////
	long nx = Inp_D->nx; //total number of cells in x direction
	long ny = Inp_D->ny; //total number of cells in y direction
	long nz = Inp_D->nz;  //total number of cells in z direction
	
	cout << "nx = " << nx << endl;
	cout << "ny = " << ny << endl;
	cout << "nz = " << nz << endl;
		
	//division of the Total Field Zone - for computing with threads
	long nr_DIV_x = Inp_D->nr_DIV_x;
	long nr_DIV_y = Inp_D->nr_DIV_y;
	long nr_DIV_z = Inp_D->nr_DIV_z;
	long nr_DIV_xORyORz = nr_DIV_x*nr_DIV_y*nr_DIV_z;
	long nr_Threads = 26 + nr_DIV_xORyORz;//total nr of threads (27 cells)
	cout << "Nr. of computing threads = " << nr_Threads << endl;

	//division for Data Save
	long Save_nr_DIV_x = Inp_D->Save_nr_DIV_x;
	long Save_nr_DIV_y = Inp_D->Save_nr_DIV_y;
	long Save_nr_DIV_z = Inp_D->Save_nr_DIV_z;
	cout << "Nr. of data saving threads = " << 3*Save_nr_DIV_x*Save_nr_DIV_y + 
		                                       3*Save_nr_DIV_y*Save_nr_DIV_z +
						                       3*Save_nr_DIV_x*Save_nr_DIV_z << endl;
	
	long nPML_x_1 = Inp_D->nPML_x_1; //[cells] the size of PML-1 zone in x
	long nPML_x_2 = Inp_D->nPML_x_2; //[cells] the size of PML-2 zone in x
	long nPML_y_1 = Inp_D->nPML_y_1; //[cells] the size of PML-1 zone in y
	long nPML_y_2 = Inp_D->nPML_y_2; //[cells] the size of PML-2 zone in y
	long nPML_z_1 = Inp_D->nPML_z_1; //[cells] the size of PML-1 zone in z
	long nPML_z_2 = Inp_D->nPML_z_2; //[cells] the size of PML-2 zone in z

	double dx = Inp_D->dx;
	double dy = Inp_D->dy;
	double dz = Inp_D->dz;
	
	//copute the time steep according to the Courant limit
	double dt = ( 1.0/speed_c/sqrt(1.0/(dx*dx)+1.0/(dy*dy)+1.0/(dz*dz)) )*0.9;
	cout << "dt = " << dt << endl;
	
	long num_iter = Inp_D->num_iter; //the number of FDTD iterations
	double limit_time = Inp_D->limit_time;//[s] 23:30 the running time, 
	                          //after it saves out all workspace data in binary

	int load_workspace_data = Inp_D->load_workspace_data; //1 - load the workspace data and continue computations
	
	int jel_limit_time = 0;

	//save out slices of the field components
	int save_field = Inp_D->save_field; //to save data - Ex, Ey, Ez, Hx, Hy, Hz slices
	long nr_Save, saveFROMinst, saveTOinst;
	long n_x_a, n_x_b, n_y_a, n_y_b, n_z_a, n_z_b; //specify the saved data 
	if (save_field == 1)
	{
	    nr_Save = Inp_D->nr_Save;
		saveFROMinst = Inp_D->saveFROMinst; //from
		saveTOinst   = Inp_D->saveTOinst;//to
		
		//specify the saving area 
		n_x_a = Inp_D->n_x_a;
		n_x_b = Inp_D->n_x_b;
		n_y_a = Inp_D->n_y_a;
		n_y_b = Inp_D->n_y_b;
		n_z_a = Inp_D->n_z_a;
		n_z_b = Inp_D->n_z_b;
	}

	//average over a rectangular volume
	long nx_a_av, nx_b_av, ny_a_av, ny_b_av, nz_a_av, nz_b_av;
	long Vol_av;
	if (Inp_D->aver_field_volume == 1)
	{
		nx_a_av = Inp_D->nx_a_av; nx_b_av = Inp_D->nx_b_av; 
		ny_a_av = Inp_D->ny_a_av; ny_b_av = Inp_D->ny_b_av; 
		nz_a_av = Inp_D->nz_a_av; nz_b_av = Inp_D->nz_b_av; 

		Vol_av = (nx_b_av - nx_a_av + 1)*(ny_b_av - ny_a_av + 1)*(nz_b_av - nz_a_av + 1);
	}
	
	//fourier transform in a rectangular volume
	long nx_a_f, nx_b_f, ny_a_f, ny_b_f, nz_a_f, nz_b_f;
	double frec_1, frec_2;
	long n_frec, n_f_Points, start_fourier;
	if (Inp_D->fourier_transf_vol == 1) 
	{
		//compute the Fourier transform in this volume
		nx_a_f = Inp_D->nx_a_f; nx_b_f = Inp_D->nx_b_f; 
		ny_a_f = Inp_D->ny_a_f; ny_b_f = Inp_D->ny_b_f; 
		nz_a_f = Inp_D->nz_a_f; nz_b_f = Inp_D->nz_b_f; 
		//the total nr of points
		n_f_Points = (nx_b_f - nx_a_f + 1)*(ny_b_f - ny_a_f + 1)*(nz_b_f - nz_a_f + 1);
		//between these frequencies
		frec_1 = Inp_D->frec_1;
		frec_2 = Inp_D->frec_2;
		n_frec = Inp_D->n_frec; //nr of frequencies
		start_fourier = Inp_D->start_fourier;
	}

	///////////////////////////////////////////////////////////////////////////////////////
	//the electromagnetic energy
	///////////////////////////////////////////////////////////////////////////////////////
	#ifdef eval_W
		//division - for computing with threads
		long nr_W_DIV_x = Inp_D->nr_W_DIV_x;
		long nr_W_DIV_y = Inp_D->nr_W_DIV_y;
		long nr_W_DIV_z = Inp_D->nr_W_DIV_z;

		long nr_Threads_W = nr_W_DIV_x*nr_W_DIV_y*nr_W_DIV_z;

		//defines the volume where the energy will be computed
		long nx_W_a = Inp_D->nx_W_a,  nx_W_b = Inp_D->nx_W_b; 
		long ny_W_a = Inp_D->ny_W_a,  ny_W_b = Inp_D->ny_W_b; 
		long nz_W_a = Inp_D->nz_W_a,  nz_W_b = Inp_D->nz_W_b;

		//the size of W
		long nx_W = nx_W_b - nx_W_a;
		long ny_W = ny_W_b - ny_W_a;
		long nz_W = nz_W_b - nz_W_a;
		
		double Vol = nx_W*ny_W*nz_W*dx*dy*dz;
		double Surf_xz = nx_W*nz_W*dx*dz;

		//compute the energy between the instants
		long eval_W_FROMinst = Inp_D->eval_W_FROMinst;
		long eval_W_TOinst   = Inp_D->eval_W_TOinst;

		//the average of the energy in the specified volume
		long n_W_atl = eval_W_TOinst - eval_W_FROMinst;
		long iter_W = 0;

		//save slices from W
		int save_W = Inp_D->save_W; //to save W
		long nr_W_Save = Inp_D->nr_W_Save; //to save at
		long save_W_FROMinst = Inp_D->save_W_FROMinst; //>= eval_W_FROMinst
		long save_W_TOinst   = Inp_D->save_W_TOinst; //<= eval_W_TOinst
		//defines the volume where the energy will be saved
		//take into consideration nx_W_a, nx_W_b, ... nz_W_b
		long n_x_W_a = Inp_D->n_x_W_a,    n_x_W_b = Inp_D->n_x_W_b; 
		long n_y_W_a = Inp_D->n_y_W_a,    n_y_W_b = Inp_D->n_y_W_b; 
		long n_z_W_a = Inp_D->n_z_W_a,    n_z_W_b = Inp_D->n_z_W_b;

		//division - for saving W with threads
		long Save_nr_W_DIV_x = Inp_D->Save_nr_W_DIV_x;
		long Save_nr_W_DIV_y = Inp_D->Save_nr_W_DIV_y;
		long Save_nr_W_DIV_z = Inp_D->Save_nr_W_DIV_z;
	#endif

	///////////////////////////////////////////////////////////////////////////////////////
	//the excitation
	///////////////////////////////////////////////////////////////////////////////////////
	int source_type = Inp_D->source_type;       // 1- Gauss; 2- Sin; 3- Gauss-Sin

	int jel_plane_wave = Inp_D->jel_plane_wave; //the type of the excitation 0 - point  source
	                                            //                           1 - plane wave
	int pt_source_Ez, pt_source_Hz;//specify the polarization
	int switch_off_time;
	double const_alfa;
		
	long n_Coord;
	long **Coord_ptSource = NULL;
	double **Param_Source = NULL;
	
	double X0, tw, t0, omega, phase, teta, phi, gamma;
	long n_TS_xa, n_TS_xb, n_TS_ya, n_TS_yb,n_TS_za, n_TS_zb;
	long n_TS = 0;

	long i;

	switch (jel_plane_wave)
	{
		case 0:
			if (Inp_D->readFile_Param_Pt_Source == 1)
			{
				/////////////////////////////////////////////////???????????????
				cout << "not implemented yet" << endl;
				return 1;
			}
			else
			{
				n_Coord = 1;
				Coord_ptSource = Init_Matrix_2D<long >(n_Coord,3);
				if (!Coord_ptSource)
				{
					cout << "Memory allocation problem - **Coord_ptSource" << endl;
					return -1;
				}

				Coord_ptSource[0][0] = nx/2+2;//175;//nx/2+2;//125;
				Coord_ptSource[0][1] = ny/2+3;//835;//522;//30
				Coord_ptSource[0][2] = nz/2;

				if (source_type == 1 || source_type == 2)
				{
					Param_Source = Init_Matrix_2D<double>(n_Coord,3);
				}
				else
				{
					Param_Source = Init_Matrix_2D<double>(n_Coord,5);
				}
				if (!Param_Source)
				{
					cout << "Memory allocation problem - **Param_Source" << endl;
					return -1;
				}
				
				X0 = Inp_D->X0; //amplituide
				switch ( Inp_D->source_type)
				{
					case 1: //Gaussian source
						tw = Inp_D->tw;
						t0 = Inp_D->t0;
						break;
					case 2: //Sin source
						omega      = Inp_D->omega;
						phase      = Inp_D->phase;
						const_alfa = Inp_D->const_alfa;
						break;
					case 3://Gaussian-Sin source
						tw    = Inp_D->tw;
						t0    = Inp_D->t0;
						omega = Inp_D->omega;