?? fdtd_1d_eh_pml_loss.cpp
字號:
#include "FDTD_1D_EH_PML_LOSS.h"
#include "run_enviroment.h"
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_1D_EH_PML_LOSS::CFDTD_1D_EH_PML_LOSS(void)
{
E_1D = NULL;
//F_1 = NULL;
F_2 = NULL;
H_1D = NULL;
K_F_a = NULL;
K_F_b = NULL;
K_H_a = NULL;
K_H_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_1D/m]
num_threads = 1;
}
CFDTD_1D_EH_PML_LOSS::~CFDTD_1D_EH_PML_LOSS(void)
{
Free_Mem();
}
///////////////////////////////////////////////////////////////////////////////////////
//Initialize the nr of threads
///////////////////////////////////////////////////////////////////////////////////////
#ifdef _OPENMP
void CFDTD_1D_EH_PML_LOSS::Init_nr_THR(int nr_threads)
{
num_threads = nr_threads;
}
#endif
///////////////////////////////////////////////////////////////////////////////////////
//Initialize the class parameters
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_1D_EH_PML_LOSS::Init(long N, long n_pml, double d_t, double D, double Eps_r,
double Mu_r, double Sigma)
{
n = N;
dt = d_t;
d = D;
inv_d = 1.0/d;
n_PML = n_pml;
n_TS_a = n_PML+10;
//some usefull constants
n_PML_MIN_1 = n_PML - 1;
n_MIN_n_PML_MIN_1 = n - n_PML - 1;
n_MIN_1 = n - 1;
n_MIN_n_PML = n - n_PML;
eps_r = Eps_r;
mu_r = Mu_r;
sigma = Sigma;
beta = sqrt(eps_0*eps_r*mu_0*mu_r);
eta = sqrt( mu_0*mu_r/(eps_0*eps_r) );
E_1D = (double *)calloc(n,sizeof(double));
if (!E_1D)
{
return false;
}
/*F_1 = (double *)calloc(n_PML,sizeof(double));
if (!F_1)
{
Free_Mem();
return false;
}*/
F_2 = (double *)calloc(n_PML,sizeof(double));
if (!F_2)
{
Free_Mem();
return false;
}
K_F_a = (double *)calloc(n_PML,sizeof(double));
if (!K_F_a)
{
Free_Mem();
return false;
}
K_F_b = (double *)calloc(n_PML,sizeof(double));
if (!K_F_b)
{
Free_Mem();
return false;
}
H_1D = (double *)calloc(n_MIN_1,sizeof(double));
if (!H_1D)
{
Free_Mem();
return false;
}
K_H_a = (double *)calloc(n_PML,sizeof(double));
if (!K_H_a)
{
Free_Mem();
return false;
}
K_H_b = (double *)calloc(n_PML,sizeof(double));
if (!K_H_b)
{
Free_Mem();
return false;
}
K_E_a = ( 2.0*eps_r*eps_0 - sigma*dt )/( 2.0*eps_r*eps_0 + sigma*dt );
K_E_b = 2.0*dt/( 2.0*eps_r*eps_0 + sigma*dt )*inv_d;
K_H = dt*inv_d/(mu_0*mu_r);
return true;
}
///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Init_PML_Param()
{
long i;
//PML parameters
double ka_max = 2.0;
double exponent = 4.0;
double R_err = 1.0e-16;
double eta = sqrt( mu_0*mu_r/(eps_0*eps_r) );
double sigma_PML, sigma_max, ka_PML;
sigma_max= -(exponent+1.0)*log(R_err)/(2.0*eta*n_PML*d);
for (i = 0; i < n_PML; i++)
{
sigma_PML = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent);
ka_PML = 1.0 + (ka_max - 1.0)*pow( (n_PML-i)/((double) n_PML) ,exponent);
K_F_a[i] = (2.0*eps_0*ka_PML - sigma_PML*dt)/(2.0*eps_0*ka_PML + sigma_PML*dt);
K_F_b[i] = 2.0*eps_0/(2.0*eps_0*ka_PML + sigma_PML*dt);
sigma_PML = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent);
ka_PML = 1.0 + (ka_max - 1.0)*pow( (n_PML - i - 0.5)/n_PML ,exponent);
K_H_a[i] = (2.0*eps_0*ka_PML - sigma_PML*dt)/(2.0*eps_0*ka_PML + sigma_PML*dt);
K_H_b[i] = 2.0*eps_0*dt/( (2.0*eps_0*ka_PML + sigma_PML*dt)*d*mu_0*mu_r);
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Init_Sin_Source(double E_0, double om, double Phase)
{
E0 = E_0;
omega = om;
phase = Phase;
alfa = 1.0*omega/(2.0*pi);
Source_Type = 1;
}
///////////////////////////////////////////////////////////////////////////////////////
//Initialize a Gaussian source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_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_EH_PML_LOSS::Init_GaussSin_Source(double E_0, double om, double Phase, double t_0, double t_w)
{
E0 = E_0;
omega = om;
phase = Phase;
t0 = t_0;
tw = t_w;
Source_Type = 3;
}
///////////////////////////////////////////////////////////////////////////////////////
//Calculate E_1D field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Calc_E(double time)
{
long i, ii, ia;
double Fz_r;
/*for (i = 1; i < n_PML; i++)
{
Fz_r = F_1[i];
F_1[i] = K_F_a[i]*F_1[i] + K_F_b[i]*( H_1D[i-1] - H_1D[i]);
E_1D[i] = K_E_a*E_1D[i] + K_E_b*( F_1[i] - Fz_r );
}*/
#pragma omp parallel for shared(H_1D,K_E_a,K_E_b,E_1D,n_MIN_n_PML) private(i) num_threads(num_threads)
for (i = 1; i < n_MIN_n_PML; i++)
{
E_1D[i] = K_E_a*E_1D[i] + K_E_b*( H_1D[i-1] - H_1D[i]);
}
ia = 0;
ii = n_PML_MIN_1;
for (i = n_MIN_n_PML; i < n_MIN_1; i++)
{
Fz_r = F_2[ia];
F_2[ia] = K_F_a[ii]*F_2[ia] + K_F_b[ii]*( H_1D[i-1] - H_1D[i]);
E_1D[i] = K_E_a*E_1D[i] + K_E_b*( F_2[ia] - Fz_r );
ia++;
ii--;
}
switch (Source_Type)
{
case 1:
E_1D[0] = (1.0-exp(-alfa*time))*E0*cos(omega*time + phase);
break;
case 2:
E_1D[0] = E0*exp( - (time - t0)*(time - t0)/(tw*tw) );
break;
case 3:
E_1D[0] = E0*cos(omega*(time-t0) + phase)*exp( - (time - t0)*(time - t0)/(tw*tw) );
break;
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Calculate H_1D field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Calc_H()
{
long i, ii;
/*for (i = 0; i < n_PML; i++)
{
H_1D[i] = K_H_a[i]*H_1D[i] + K_H_b[i]*( E_1D[i] - E_1D[i+1]);
}*/
#pragma omp parallel for shared(H_1D,K_H,E_1D,n_MIN_n_PML_MIN_1) private(i) num_threads(num_threads)
for (i = 0; i < n_MIN_n_PML_MIN_1; i++)
{
H_1D[i] += K_H*( E_1D[i] - E_1D[i+1]);
}
ii = n_PML_MIN_1;
for (i = n_MIN_n_PML_MIN_1; i < n_MIN_1; i++)
{
H_1D[i] = K_H_a[ii]*H_1D[i] + K_H_b[ii]*( E_1D[i] - E_1D[i+1] );
ii--;
}
}
///////////////////////////////////////////////////////////////////////////////////////
//Get the address of E_1D and H_1D
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Get_Data(double *&e, double *&h)
{
e = E_1D;
h = H_1D;
}
///////////////////////////////////////////////////////////////////////////////////////
//Save out the workspace data - the possibility to continue computations
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_1D_EH_PML_LOSS::Save_FDTD_1D_Workspace(char *path)
{
char *path_name = NULL, *file_name = NULL;
path_name =(char *) calloc(256,sizeof(char));
file_name =(char *) calloc(256,sizeof(char));
if (!path_name && !file_name)
{
if (path_name)
free(path_name);
if (file_name)
free(file_name);
return 1;
}
strcpy(path_name,path);
strcat(path_name,"/data_fdtd_1D");
//make directory to save F G E_1D
#ifdef run_enviroment
CreateDirectory(path_name,0);
#else
mode_t Mode = S_IRWXU;
mkdir(path_name,Mode);//for UNIX-AIX
#endif
//save F_2
strcpy(file_name,path_name);
strcat(file_name,"/F_1D");
if (save_1D_binary(F_2, n_PML, 0, file_name))
{
free(path_name);
free(file_name);
return 2; //faild to save data
}
//save E_1D
strcpy(file_name,path_name);
strcat(file_name,"/E_1D");
if (save_1D_binary(E_1D, n, 0, file_name))
{
free(path_name);
free(file_name);
return 2; //faild to save data
}
//save H_1D
strcpy(file_name,path_name);
strcat(file_name,"/H_1D");
if(save_1D_binary(H_1D, n-1, 0, file_name))
{
free(path_name);
free(file_name);
return 2; //faild to save data
}
free(path_name);
free(file_name);
return 0;
}
int CFDTD_1D_EH_PML_LOSS::Load_FDTD_1D_Workspace(char *path)
{
char *file_name = NULL;
file_name =(char *) calloc(256,sizeof(char));
if (!file_name)
{
return 1;
}
//load F2
strcpy(file_name,path);
strcat(file_name,"/F_1D_0.dat");
if (load_1D_binary(F_2, n_PML, file_name))
{
free(file_name);
return 2; //faild to load data
}
//load E_1D
strcpy(file_name,path);
strcat(file_name,"/E_1D_0.dat");
if (load_1D_binary(E_1D, n, file_name))
{
free(file_name);
return 2; //faild to load data
}
//load Hx_1D
strcpy(file_name,path);
strcat(file_name,"/H_1D_0.dat");
if (load_1D_binary(H_1D, n-1, file_name))
{
free(file_name);
return 2; //faild to load data
}
free(file_name);
return 0;
}
///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_PML_LOSS::Free_Mem(void)
{
if (E_1D)
{
free(E_1D);
E_1D = NULL;
}
/*if (F_1)
{
free(F_1);
F_1 = NULL;
}*/
if (F_2)
{
free(F_2);
F_2 = NULL;
}
if (K_F_a)
{
free(K_F_a);
K_F_a = NULL;
}
if (K_F_b)
{
free(K_F_b);
K_F_b = NULL;
}
if (H_1D)
{
free(H_1D);
H_1D = NULL;
}
if (K_H_a)
{
free(K_H_a);
K_H_a = NULL;
}
if (K_H_b)
{
free(K_H_b);
K_H_b = NULL;
}
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -