fdtd_1d_ezhy_pml_loss.cpp

2维时域有限差分模拟

#include "StdAfx.h"
#include ".\fdtd_1d_ezhy_pml_loss.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_1D_EZHY_PML_LOSS::CFDTD_1D_EZHY_PML_LOSS(void)
{
	Ez = NULL;
	Fz = NULL;
	Hy = NULL;

	K_Fz_a = NULL;
	K_Fz_b = NULL;
	K_Hy_a = NULL;
	K_Hy_b = NULL;

	pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
	//permittivity of free space 
	eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]
}

CFDTD_1D_EZHY_PML_LOSS::~CFDTD_1D_EZHY_PML_LOSS(void)
{
	Free_Mem();
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the class parameters
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_1D_EZHY_PML_LOSS::Init(int n_x, int n_pml, double d_t, double d_x, double Eps_r, 
								  double Mu_r, double Sigma)
{
	nx = n_x;
	dt = d_t;
	dx = d_x;
	n_PML = n_pml;

	eps_r = Eps_r;
	mu_r = Mu_r;
	sigma = Sigma;


	Ez = (double *)calloc(nx,sizeof(double));
	if (!Ez)
	{
		return false;
	}

	Fz = (double *)calloc(nx,sizeof(double));
	if (!Fz)
	{
		Free_Mem();
		return false;
	}

	K_Fz_a = (double *)calloc(nx,sizeof(double));
	if (!K_Fz_a)
	{
		Free_Mem();
		return false;
	}

	K_Fz_b = (double *)calloc(nx,sizeof(double));
	if (!K_Fz_b)
	{
		Free_Mem();
		return false;
	}

	Hy = (double *)calloc(nx-1,sizeof(double));
	if (!Hy)
	{
		Free_Mem();
		return false;
	}

	K_Hy_a = (double *)calloc(nx-1,sizeof(double));
	if (!K_Hy_a)
	{
		Free_Mem();
		return false;
	}

 	K_Hy_b = (double *)calloc(nx-1,sizeof(double));
	if (!K_Hy_b)
	{
		Free_Mem();
		return false;
	}

	K_Ez_a = ( 2.0*eps_r*eps_0 - sigma*dt )/( 2.0*eps_r*eps_0 + sigma*dt );
	K_Ez_b =  2.0/( 2.0*eps_r*eps_0 + sigma*dt );
	
	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::Init_PML_Param(void)
{
	int i;

	for (i = 0; i<nx-1; i++)
	{
		K_Fz_a[i] = 1.0;
		K_Fz_b[i] = dt/dx;

		K_Hy_a[i] = 1.0;
		K_Hy_b[i] = dt/(mu_0*mu_r*dx);
	}
	K_Fz_a[nx-1] = 1.0;
	K_Fz_b[nx-1] = dt/dx;

	//PML parameters
	double ka_max = 2.0;
	int exponent = 4;
	double R_err = 1.0e-16;
	double eta = sqrt( mu_0*mu_r/(eps_0*eps_r) );

	double sigma_x, sigma_max, ka_x;
		
	sigma_max= -(exponent+1.0)*log(R_err)/(2.0*eta*n_PML*dx);

	for (i = 0; i<n_PML; i++)
	{
		sigma_x         = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent);
		ka_x            = 1.0 + (ka_max - 1.0)*pow( (n_PML-i)/((double) n_PML) ,exponent);

		K_Fz_a[i]       = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*ka_x + sigma_x*dt);
		K_Fz_a[nx-i-1]  = K_Fz_a[i];

		K_Fz_b[i]       = 2.0*eps_0*dt/(2.0*eps_0*ka_x + sigma_x*dt)/dx;
		K_Fz_b[nx-i-1]  = K_Fz_b[i];

		sigma_x         = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		ka_x            = 1.0 + (ka_max - 1.0)*pow( (n_PML - i - 0.5)/n_PML ,exponent);

        K_Hy_a[i]       = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*ka_x + sigma_x*dt);
		K_Hy_a[nx-i-2]  = K_Hy_a[i];

		K_Hy_b[i]       = 2.0*eps_0*dt/( (2.0*eps_0*ka_x + sigma_x*dt)*dx*mu_0*mu_r);
		K_Hy_b[nx-i-2]  = K_Hy_b[i];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::Init_Sin_Source(double E_0, double om, double Phi)
{
	E0 = E_0;
	omega = om;
	phi = Phi;

	alfa = 3*omega/(2*pi);

	Source_Type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize  a Gaussian source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::Init_Gauss_Source(double E_0, double t_0, double t_w)
{
	E0 = E_0;
	t0 = t_0;
	tw = t_w;
	
	Source_Type = 2;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal modulated Gaussian pulse
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::Init_GaussSin(double E_0, double om, double Phi, double t_0, double t_w)
{
	E0 = E_0;
	omega = om;
	phi = Phi;
	t0 = t_0;
	tw = t_w;

	Source_Type = 3;
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ez field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::calc_Ez(double time)
{
	int i;
	double Fz_r;

	for (i = 1; i < nx-1; i++)
	{
		Fz_r = Fz[i];
		
		Fz[i] = K_Fz_a[i]*Fz[i] + K_Fz_b[i]*( Hy[i] - Hy[i-1] );

		Ez[i] = K_Ez_a*Ez[i] + K_Ez_b*( Fz[i] - Fz_r );
	}

	switch (Source_Type)
	{
		case 1:
			Ez[n_PML+1] = (1-exp(-alfa*time))*E0*cos(omega*time + phi);
			break;
		case 2:
			Ez[n_PML+1] = E0*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
		case 3:
			Ez[n_PML+1] = E0*cos(omega*(time-t0) + phi)*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hy field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::calc_Hy()
{
	int i;

	for (i = 0; i < nx-1; i++)
	{
		Hy[i] = K_Hy_a[i]*Hy[i] + K_Hy_b[i]*( Ez[i+1] - Ez[i] );
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Get the address of Ez and Hy
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::Get_Data(double *&ez, double *&hy)
{
	ez = Ez;
	hy = Hy;
}

///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EZHY_PML_LOSS::Free_Mem(void)
{
	if (Ez)
	{
		free(Ez);
		Ez = NULL;
	}

	if (Fz)
	{
		free(Fz);
		Fz = NULL;
	}

	if (K_Fz_a)
	{
		free(K_Fz_a);
		K_Fz_a = NULL;
	}

	if (K_Fz_b)
	{
		free(K_Fz_b);
		K_Fz_b = NULL;
	}

	if (Hy)
	{
		free(Hy);
		Hy = NULL;
	}

	if (K_Hy_a)
	{
		free(K_Hy_a);
		K_Hy_a = NULL;
	}

	if (K_Hy_b)
	{
		free(K_Hy_b);
		K_Hy_b = NULL;
	}
}