diff --git a/PlotScript/Animations/animation_Bx.gif b/PlotScript/Animations/animation_Bx.gif deleted file mode 100644 index dc08b37..0000000 Binary files a/PlotScript/Animations/animation_Bx.gif and /dev/null differ diff --git a/PlotScript/Animations/animation_By.gif b/PlotScript/Animations/animation_By.gif deleted file mode 100644 index 5bc2416..0000000 Binary files a/PlotScript/Animations/animation_By.gif and /dev/null differ diff --git a/PlotScript/Animations/animation_Bz.gif b/PlotScript/Animations/animation_Bz.gif deleted file mode 100644 index 31d94b6..0000000 Binary files a/PlotScript/Animations/animation_Bz.gif and /dev/null differ diff --git a/PlotScript/Animations/animation_Ey.gif b/PlotScript/Animations/animation_Ey.gif deleted file mode 100644 index fb33b3d..0000000 Binary files a/PlotScript/Animations/animation_Ey.gif and /dev/null differ diff --git a/PlotScript/Animations/animation_Ez.gif b/PlotScript/Animations/animation_Ez.gif deleted file mode 100644 index 1cfb852..0000000 Binary files a/PlotScript/Animations/animation_Ez.gif and /dev/null differ diff --git a/PlotScript/src/Release/sample.exe b/PlotScript/src/Release/sample.exe index bdef30a..d23c400 100644 Binary files a/PlotScript/src/Release/sample.exe and b/PlotScript/src/Release/sample.exe differ diff --git a/include/FDTD.h b/include/FDTD.h index d70e355..ead71e2 100644 --- a/include/FDTD.h +++ b/include/FDTD.h @@ -15,7 +15,6 @@ class Field private: int Ni, Nj, Nk; std::vector field; - double nul = 0.0; public: Field(const int _Ni, const int _Nj, const int _Nk); @@ -33,36 +32,34 @@ class FDTD private: Field Ex, Ey, Ez, Bx, By, Bz; - Field Exy, Exz, Eyx, Eyz, Ezx, Ezy; - Field Bxy, Bxz, Byx, Byz, Bzx, Bzy; - Field EsigmaX, EsigmaY, EsigmaZ, BsigmaX, BsigmaY, BsigmaZ; - std::vector Jx; std::vector Jy; std::vector Jz; - /*std::vector EsigmaX; - std::vector EsigmaY; - std::vector EsigmaZ; - std::vector BsigmaX; - std::vector BsigmaY; - std::vector BsigmaZ;*/ + Field Exy, Exz, Eyx, Eyz, Ezx, Ezy; + Field Bxy, Bxz, Byx, Byz, Bzx, Bzy; + Field EsigmaX, EsigmaY, EsigmaZ, + BsigmaX, BsigmaY, BsigmaZ; Parameters parameters; - PMLparameters pml; - double pml_percent; + double pml_percent, dt; int pml_size_i, pml_size_j, pml_size_k; - double dt; void update_E(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t); - void update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t); + void update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2]); - double PMLcoef(double sigma, double constant); + double PMLcoef(double sigma); void update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]); void update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]); - void set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], double SGm[2], std::function dist); - void set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2], double SGm[2], std::function dist); - void set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], double SGm[2], std::function dist); + + void set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], + double SGm, std::function dist); + void set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2], + double SGm, std::function dist); + void set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], + double SGm, std::function dist); + + double calc_max_value(Field& field); public: FDTD(Parameters _parameters, double _dt, double _pml_percent); @@ -71,4 +68,7 @@ class FDTD std::vector& get_current(Component); std::vector> update_fields(const int time); + + double calc_reflection(Field E_start[3], Field B_start[3], + Field E_final[3], Field B_final[3]); }; diff --git a/include/Structures.h b/include/Structures.h index 359e448..28f21d4 100644 --- a/include/Structures.h +++ b/include/Structures.h @@ -5,17 +5,15 @@ namespace FDTDconst { const double C = 3e10; - const double R = 0.01; - const double EPS0 = 1.60217656535e-12; - const double PI = 3.1415926535; - const double MU0 = 1e6 * EPS0; + const double R = 1e-12; + const double EPS0 = 1.0; + const double MU0 = EPS0; const double N = 4.0; } namespace FDTDstruct { enum class Component { EX, EY, EZ, BX, BY, BZ, JX, JY, JZ }; - enum class Axis { X, Y, Z }; struct SelectedFields { @@ -33,15 +31,6 @@ namespace FDTDstruct double period_z = static_cast(m) * FDTDconst::C; }; - struct PMLparameters { - double sigma_low_xy; - double sigma_low_yz; - double sigma_low_zx; - double sigma_up_xy; - double sigma_up_yz; - double sigma_up_zx; - }; - struct Parameters { int Ni; int Nj; @@ -49,10 +38,8 @@ namespace FDTDstruct double ax; double bx; - double ay; double by; - double az; double bz; diff --git a/samples/sample.cpp b/samples/sample.cpp index eb7a499..0ea08cf 100644 --- a/samples/sample.cpp +++ b/samples/sample.cpp @@ -36,7 +36,8 @@ Axis get_axis(Component field_E, Component field_B) return selected_axis; } -void plane_wave(std::vector arguments, std::vector numbers, char* outfile_path) +void plane_wave(std::vector arguments, + std::vector numbers, char* outfile_path) { // Saving data from the console Component selected_field_1 = static_cast(std::atoi(arguments[1])); @@ -53,17 +54,13 @@ void plane_wave(std::vector arguments, std::vector numbers, char* double sizes_z[2] = { numbers[7], numbers[8] }; // { az, bz } int iterations_num = static_cast(numbers[9]); double time = numbers[10]; - double time_step = 0.4; // time / static_cast(iterations_num); // dt + double time_step = time / static_cast(iterations_num); // dt // Initialization of the initializing function and the true solution function std::function initial_func = [](double x, double size[2]) { return sin(2.0 * M_PI * (x - size[0]) / (size[1] - size[0])); }; - std::function true_func = [](double x, double t, double size[2]) - { - return sin(2.0 * M_PI * (x - size[0] - FDTDconst::C * t) / (size[1] - size[0])); - }; // Determination of the wave propagation axis Axis axis = get_axis(flds_selected[0], flds_selected[1]); @@ -80,17 +77,17 @@ void plane_wave(std::vector arguments, std::vector numbers, char* grid_sizes[0], grid_sizes[1], grid_sizes[2], - 0.0, + sizes_x[0], sizes_x[1], - 0.0, + sizes_y[0], sizes_y[1], - 0.0, + sizes_z[0], sizes_z[1], - FDTDconst::C, //(sizes_x[1] - sizes_x[0]) / static_cast(grid_sizes[0]), - FDTDconst::C, //(sizes_y[1] - sizes_y[0]) / static_cast(grid_sizes[1]), - FDTDconst::C, //(sizes_z[1] - sizes_z[0]) / static_cast(grid_sizes[2]) + (sizes_x[1] - sizes_x[0]) / static_cast(grid_sizes[0]), + (sizes_y[1] - sizes_y[0]) / static_cast(grid_sizes[1]), + (sizes_z[1] - sizes_z[0]) / static_cast(grid_sizes[2]) }; - FDTD method(params, time_step, 0.1); + FDTD method(params, time_step, 0.0); // Meaningful calculations Test_FDTD test(params); @@ -104,15 +101,13 @@ void plane_wave(std::vector arguments, std::vector numbers, char* } std::vector> data = method.update_fields(iterations_num); - - //std::cout << test.get_max_abs_error(method.get_field(fld_tested), fld_tested, true_func, time) << std::endl; if (write_flag == true) { // Writing the results to a file try { - write_spherical(data, Axis::X, iterations_num, outfile_path); + write_plane(method, axis, outfile_path); } catch (const std::exception& except) { @@ -125,18 +120,22 @@ void spherical_wave(int n, int it, char* base_path = "") { CurrentParameters cur_param { - 2, + 8, 4, - 0.4, + 0.2, }; double T = cur_param.period; double Tx = cur_param.period_x; double Ty = cur_param.period_y; double Tz = cur_param.period_z; cur_param.iterations = static_cast(static_cast(cur_param.period) / cur_param.dt); - std::function cur_func = [T, Tx, Ty, Tz](double x, double y, double z, double t) + std::function cur_func + = [T, Tx, Ty, Tz](double x, double y, double z, double t) { - return sin(2.0 * M_PI * t / T) * pow(cos(2.0 * M_PI * x / Tx), 2.0) * pow(cos(2.0 * M_PI * y / Ty), 2.0) * pow(cos(2.0 * M_PI * z / Tz), 2.0); + return sin(2.0 * M_PI * t / T) + * pow(cos(2.0 * M_PI * x / Tx), 2.0) + * pow(cos(2.0 * M_PI * y / Ty), 2.0) + * pow(cos(2.0 * M_PI * z / Tz), 2.0); }; // Initialization of the structures and method @@ -182,8 +181,8 @@ int main(int argc, char* argv[]) if (argc == 1) { #ifdef __USE_SPHERICAL_WAVE__ - int N = 50; // 100; - int Iterations = 250; // 200; + int N = 70; + int Iterations = 300; spherical_wave(N, Iterations, "../../PlotScript/"); #endif diff --git a/src/FDTD.cpp b/src/FDTD.cpp index eea3e8d..8f493d5 100644 --- a/src/FDTD.cpp +++ b/src/FDTD.cpp @@ -1,7 +1,12 @@ #include "FDTD.h" + #include -Field::Field(const int _Ni = 1, const int _Nj = 1, const int _Nk = 1) : Ni(_Ni), Nj(_Nj), Nk(_Nk) +#define __DEBUG_PML__ +//#define __OUTPUT_SIGMA__ + +Field::Field(const int _Ni = 1, const int _Nj = 1, const int _Nk = 1) + : Ni(_Ni), Nj(_Nj), Nk(_Nk) { int size = Ni * Nj * Nk; field = std::vector(size, 0.0); @@ -28,16 +33,19 @@ double& Field::operator() (int i, int j, int k) int j_isNj = (j == Nj); int k_isNk = (k == Nk); - int truly_i = (Ni - 1) * i_isMinusOne + i * !(i_isMinusOne || i_isNi); - int truly_j = (Nj - 1) * j_isMinusOne + j * !(j_isMinusOne || j_isNj); - int truly_k = (Nk - 1) * k_isMinusOne + k * !(k_isMinusOne || k_isNk); + int truly_i = (Ni - 1) * i_isMinusOne + i * + !(i_isMinusOne || i_isNi); + int truly_j = (Nj - 1) * j_isMinusOne + j * + !(j_isMinusOne || j_isNj); + int truly_k = (Nk - 1) * k_isMinusOne + k * + !(k_isMinusOne || k_isNk); int index = truly_i + truly_j * Ni + truly_k * Ni * Nj; return field[index]; } void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], - double SGm[2], std::function dist) + double SGm, std::function dist) { #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) @@ -46,16 +54,18 @@ void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - EsigmaX(i, j, k) = SGm[0] * pow((static_cast(dist(i, j, k))) - / static_cast(pml_size_i), FDTDconst::N); - BsigmaX(i, j, k) = SGm[1] * pow((static_cast(dist(i, j, k)) + 0.5) - / static_cast(pml_size_i), FDTDconst::N); + EsigmaX(i, j, k) = SGm * + pow((static_cast(dist(i, j, k))) / + static_cast(pml_size_i), FDTDconst::N); + BsigmaX(i, j, k) = SGm * + pow((static_cast(dist(i, j, k)) + 0.5) / + static_cast(pml_size_i), FDTDconst::N); } } } } void FDTD::set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2], - double SGm[2], std::function dist) + double SGm, std::function dist) { #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) @@ -64,16 +74,18 @@ void FDTD::set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2], { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - EsigmaY(i, j, k) = SGm[0] * pow((static_cast(dist(i, j, k))) - / static_cast(pml_size_j), FDTDconst::N); - BsigmaY(i, j, k) = SGm[1] * pow((static_cast(dist(i, j, k)) + 0.5) - / static_cast(pml_size_j), FDTDconst::N); + EsigmaY(i, j, k) = SGm * + pow((static_cast(dist(i, j, k))) / + static_cast(pml_size_j), FDTDconst::N); + BsigmaY(i, j, k) = SGm * + pow((static_cast(dist(i, j, k)) + 0.5) / + static_cast(pml_size_j), FDTDconst::N); } } } } void FDTD::set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], - double SGm[2], std::function dist) + double SGm, std::function dist) { #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) @@ -82,17 +94,17 @@ void FDTD::set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - EsigmaZ(i, j, k) = SGm[0] * pow((static_cast(dist(i, j, k))) + EsigmaZ(i, j, k) = SGm * + pow((static_cast(dist(i, j, k))) / static_cast(pml_size_k), FDTDconst::N); - BsigmaZ(i, j, k) = SGm[1] * pow((static_cast(dist(i, j, k)) + 0.5) + BsigmaZ(i, j, k) = SGm * + pow((static_cast(dist(i, j, k)) + 0.5) / static_cast(pml_size_k), FDTDconst::N); } } } } - - FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) : parameters(_parameters), dt(_dt), pml_percent(_pml_percent) { @@ -109,10 +121,14 @@ FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) : { throw std::exception("ERROR: invalid parameters"); } - Ex = Ey = Ez = Bx = By = Bz = Field(parameters.Ni, parameters.Nj, parameters.Nk); - Exy = Exz = Eyx = Eyz = Ezx = Ezy = Field(parameters.Ni, parameters.Nj, parameters.Nk); - Bxy = Bxz = Byx = Byz = Bzx = Bzy = Field(parameters.Ni, parameters.Nj, parameters.Nk); - EsigmaX = EsigmaY = EsigmaZ = BsigmaX = BsigmaY = BsigmaZ = Field(parameters.Ni, parameters.Nj, parameters.Nk); + Ex = Ey = Ez = Bx = By = Bz + = Field(parameters.Ni, parameters.Nj, parameters.Nk); + Exy = Exz = Eyx = Eyz = Ezx = Ezy + = Field(parameters.Ni, parameters.Nj, parameters.Nk); + Bxy = Bxz = Byx = Byz = Bzx = Bzy + = Field(parameters.Ni, parameters.Nj, parameters.Nk); + EsigmaX = EsigmaY = EsigmaZ = BsigmaX = BsigmaY = BsigmaZ + = Field(parameters.Ni, parameters.Nj, parameters.Nk); pml_size_i = static_cast(static_cast(parameters.Ni) * pml_percent); pml_size_j = static_cast(static_cast(parameters.Nj) * pml_percent); @@ -166,15 +182,21 @@ void FDTD::update_E(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t) { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - Ex(i, j, k) = Ex(i, j, k) - 4.0 * M_PI * dt * Jx[t](i, j, k) + FDTDconst::C * dt * ((Bz(i, j, k) - Bz(i, j - 1, k)) / dy - (By(i, j, k) - By(i, j, k - 1)) / dz); - Ey(i, j, k) = Ey(i, j, k) - 4.0 * M_PI * dt * Jy[t](i, j, k) + FDTDconst::C * dt * ((Bx(i, j, k) - Bx(i, j, k - 1)) / dz - (Bz(i, j, k) - Bz(i - 1, j, k)) / dx); - Ez(i, j, k) = Ez(i, j, k) - 4.0 * M_PI * dt * Jz[t](i, j, k) + FDTDconst::C * dt * ((By(i, j, k) - By(i - 1, j, k)) / dx - (Bx(i, j, k) - Bx(i, j - 1, k)) / dy); + Ex(i, j, k) = Ex(i, j, k) - 4.0 * M_PI * dt * Jx[t](i, j, k) + + FDTDconst::C * dt * ((Bz(i, j, k) - Bz(i, j - 1, k)) / dy - + (By(i, j, k) - By(i, j, k - 1)) / dz); + Ey(i, j, k) = Ey(i, j, k) - 4.0 * M_PI * dt * Jy[t](i, j, k) + + FDTDconst::C * dt * ((Bx(i, j, k) - Bx(i, j, k - 1)) / dz - + (Bz(i, j, k) - Bz(i - 1, j, k)) / dx); + Ez(i, j, k) = Ez(i, j, k) - 4.0 * M_PI * dt * Jz[t](i, j, k) + + FDTDconst::C * dt * ((By(i, j, k) - By(i - 1, j, k)) / dx - + (Bx(i, j, k) - Bx(i, j - 1, k)) / dy); } } } } -void FDTD::update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t) +void FDTD::update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2]) { double dx = parameters.dx; double dy = parameters.dy; @@ -187,17 +209,23 @@ void FDTD::update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t) { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - Bx(i, j, k) += FDTDconst::C * dt / 2.0 * ((Ey(i, j, k + 1) - Ey(i, j, k)) / dz - (Ez(i, j + 1, k) - Ez(i, j, k)) / dy); - By(i, j, k) += FDTDconst::C * dt / 2.0 * ((Ez(i + 1, j, k) - Ez(i, j, k)) / dx - (Ex(i, j, k + 1) - Ex(i, j, k)) / dz); - Bz(i, j, k) += FDTDconst::C * dt / 2.0 * ((Ex(i, j + 1, k) - Ex(i, j, k)) / dy - (Ey(i + 1, j, k) - Ey(i, j, k)) / dx); + Bx(i, j, k) += FDTDconst::C * dt / 2.0 * + ((Ey(i, j, k + 1) - Ey(i, j, k)) / dz - + (Ez(i, j + 1, k) - Ez(i, j, k)) / dy); + By(i, j, k) += FDTDconst::C * dt / 2.0 * + ((Ez(i + 1, j, k) - Ez(i, j, k)) / dx - + (Ex(i, j, k + 1) - Ex(i, j, k)) / dz); + Bz(i, j, k) += FDTDconst::C * dt / 2.0 * + ((Ex(i, j + 1, k) - Ex(i, j, k)) / dy - + (Ey(i + 1, j, k) - Ey(i, j, k)) / dx); } } } } -double FDTD::PMLcoef(double sigma, double constant) +double FDTD::PMLcoef(double sigma) { - return exp(-sigma * dt * FDTDconst::C);// / constant); + return exp(-sigma * dt * FDTDconst::C); } void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) @@ -205,6 +233,7 @@ void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) double dx = parameters.dx; double dy = parameters.dy; double dz = parameters.dz; + #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) { @@ -214,29 +243,29 @@ void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) { if (EsigmaX(i, j, k) != 0) { - Eyx(i, j, k) = Eyx(i, j, k) * PMLcoef(EsigmaX(i, j, k), FDTDconst::EPS0) - - (1.0 - PMLcoef(EsigmaX(i, j, k), FDTDconst::EPS0)) / (EsigmaX(i, j, k) * dx) * //FDTDconst::C * + Eyx(i, j, k) = Eyx(i, j, k) * PMLcoef(EsigmaX(i, j, k)) - + (1.0 - PMLcoef(EsigmaX(i, j, k))) / (EsigmaX(i, j, k) * dx) * (Bz(i, j, k) - Bz(i - 1, j, k)); - Ezx(i, j, k) = Ezx(i, j, k) * PMLcoef(EsigmaX(i, j, k), FDTDconst::EPS0) + - (1.0 - PMLcoef(EsigmaX(i, j, k), FDTDconst::EPS0)) / (EsigmaX(i, j, k) * dx) * //FDTDconst::C * + Ezx(i, j, k) = Ezx(i, j, k) * PMLcoef(EsigmaX(i, j, k)) + + (1.0 - PMLcoef(EsigmaX(i, j, k))) / (EsigmaX(i, j, k) * dx) * (By(i, j, k) - By(i - 1, j, k)); } if (EsigmaY(i, j, k) != 0) { - Exy(i, j, k) = Exy(i, j, k) * PMLcoef(EsigmaY(i, j, k), FDTDconst::EPS0) + - (1.0 - PMLcoef(EsigmaY(i, j, k), FDTDconst::EPS0)) / (EsigmaY(i, j, k) * dy) * //FDTDconst::C * + Exy(i, j, k) = Exy(i, j, k) * PMLcoef(EsigmaY(i, j, k)) + + (1.0 - PMLcoef(EsigmaY(i, j, k))) / (EsigmaY(i, j, k) * dy) * (Bz(i, j, k) - Bz(i, j - 1, k)); - Ezy(i, j, k) = Ezy(i, j, k) * PMLcoef(EsigmaY(i, j, k), FDTDconst::EPS0) - - (1.0 - PMLcoef(EsigmaY(i, j, k), FDTDconst::EPS0)) / (EsigmaY(i, j, k) * dy) * //FDTDconst::C * + Ezy(i, j, k) = Ezy(i, j, k) * PMLcoef(EsigmaY(i, j, k)) - + (1.0 - PMLcoef(EsigmaY(i, j, k))) / (EsigmaY(i, j, k) * dy) * (Bx(i, j, k) - Bx(i, j - 1, k)); } if (EsigmaZ(i, j, k) != 0) { - Exz(i, j, k) = Exz(i, j, k) * PMLcoef(EsigmaZ(i, j, k), FDTDconst::EPS0) - - (1.0 - PMLcoef(EsigmaZ(i, j, k), FDTDconst::EPS0)) / (EsigmaZ(i, j, k) * dz) * //FDTDconst::C * + Exz(i, j, k) = Exz(i, j, k) * PMLcoef(EsigmaZ(i, j, k)) - + (1.0 - PMLcoef(EsigmaZ(i, j, k))) / (EsigmaZ(i, j, k) * dz) * (By(i, j, k) - By(i, j, k - 1)); - Eyz(i, j, k) = Eyz(i, j, k) * PMLcoef(EsigmaZ(i, j, k), FDTDconst::EPS0) + - (1.0 - PMLcoef(EsigmaZ(i, j, k), FDTDconst::EPS0)) / (EsigmaZ(i, j, k) * dz) * //FDTDconst::C * + Eyz(i, j, k) = Eyz(i, j, k) * PMLcoef(EsigmaZ(i, j, k)) + + (1.0 - PMLcoef(EsigmaZ(i, j, k))) / (EsigmaZ(i, j, k) * dz) * (Bx(i, j, k) - Bx(i, j, k - 1)); } Ex(i, j, k) = Exz(i, j, k) + Exy(i, j, k); @@ -262,29 +291,29 @@ void FDTD::update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) { if (BsigmaX(i, j, k) != 0.0) { - Byx(i, j, k) = Byx(i, j, k) * PMLcoef(BsigmaX(i, j, k), FDTDconst::MU0) + - (1.0 - PMLcoef(BsigmaX(i, j, k), FDTDconst::MU0)) / (BsigmaX(i, j, k) * dx) * //FDTDconst::C * + Byx(i, j, k) = Byx(i, j, k) * PMLcoef(BsigmaX(i, j, k)) + + (1.0 - PMLcoef(BsigmaX(i, j, k))) / (BsigmaX(i, j, k) * dx) * (Ez(i + 1, j, k) - Ez(i, j, k)); - Bzx(i, j, k) = Bzx(i, j, k) * PMLcoef(BsigmaX(i, j, k), FDTDconst::MU0) - - (1.0 - PMLcoef(BsigmaX(i, j, k), FDTDconst::MU0)) / (BsigmaX(i, j, k) * dx) * //FDTDconst::C * + Bzx(i, j, k) = Bzx(i, j, k) * PMLcoef(BsigmaX(i, j, k)) - + (1.0 - PMLcoef(BsigmaX(i, j, k))) / (BsigmaX(i, j, k) * dx) * (Ey(i + 1, j, k) - Ey(i, j, k)); } if (BsigmaY(i, j, k) != 0.0) { - Bxy(i, j, k) = Bxy(i, j, k) * PMLcoef(BsigmaY(i, j, k), FDTDconst::MU0) - - (1.0 - PMLcoef(BsigmaY(i, j, k), FDTDconst::MU0)) / (BsigmaY(i, j, k) * dy) * //FDTDconst::C * + Bxy(i, j, k) = Bxy(i, j, k) * PMLcoef(BsigmaY(i, j, k)) - + (1.0 - PMLcoef(BsigmaY(i, j, k))) / (BsigmaY(i, j, k) * dy) * (Ez(i, j + 1, k) - Ez(i, j, k)); - Bzy(i, j, k) = Bzy(i, j, k) * PMLcoef(BsigmaY(i, j, k), FDTDconst::MU0) + - (1.0 - PMLcoef(BsigmaY(i, j, k), FDTDconst::MU0)) / (BsigmaY(i, j, k) * dy) * //FDTDconst::C * + Bzy(i, j, k) = Bzy(i, j, k) * PMLcoef(BsigmaY(i, j, k)) + + (1.0 - PMLcoef(BsigmaY(i, j, k))) / (BsigmaY(i, j, k) * dy) * (Ex(i, j + 1, k) - Ex(i, j, k)); } if (BsigmaZ(i, j, k) != 0.0) { - Bxz(i, j, k) = Bxz(i, j, k) * PMLcoef(BsigmaZ(i, j, k), FDTDconst::MU0) + - (1.0 - PMLcoef(BsigmaZ(i, j, k), FDTDconst::MU0)) / (BsigmaZ(i, j, k) * dz) * //FDTDconst::C * + Bxz(i, j, k) = Bxz(i, j, k) * PMLcoef(BsigmaZ(i, j, k)) + + (1.0 - PMLcoef(BsigmaZ(i, j, k))) / (BsigmaZ(i, j, k) * dz) * (Ey(i, j, k + 1) - Ey(i, j, k)); - Byz(i, j, k) = Byz(i, j, k) * PMLcoef(BsigmaZ(i, j, k), FDTDconst::MU0) - - (1.0 - PMLcoef(BsigmaZ(i, j, k), FDTDconst::MU0)) / (BsigmaZ(i, j, k) * dz) * //FDTDconst::C * + Byz(i, j, k) = Byz(i, j, k) * PMLcoef(BsigmaZ(i, j, k)) - + (1.0 - PMLcoef(BsigmaZ(i, j, k))) / (BsigmaZ(i, j, k) * dz) * (Ex(i, j, k + 1) - Ex(i, j, k)); } Bx(i, j, k) = Bxy(i, j, k) + Bxz(i, j, k); @@ -295,6 +324,60 @@ void FDTD::update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) } } +double FDTD::calc_reflection(Field E_start[3], Field B_start[3], + Field E_final[3], Field B_final[3]) +{ + double energy_start = 0.0; + double energy_final = 0.0; + +#pragma omp parallel for collapse(2) + for (int i = pml_size_i; i < parameters.Ni - pml_size_i; i++) + { + for (int j = pml_size_j; j < parameters.Nj - pml_size_j; j++) + { + for (int k = pml_size_k; k < parameters.Nk - pml_size_k; k++) + { + energy_start += pow(E_start[0](i, j, k), 2.0) + + pow(E_start[1](i, j, k), 2.0) + + pow(E_start[2](i, j, k), 2.0) + + pow(B_start[0](i, j, k), 2.0) + + pow(B_start[1](i, j, k), 2.0) + + pow(B_start[2](i, j, k), 2.0); + + energy_final += pow(E_final[0](i, j, k), 2.0) + + pow(E_final[1](i, j, k), 2.0) + + pow(E_final[2](i, j, k), 2.0) + + pow(B_final[0](i, j, k), 2.0) + + pow(B_final[1](i, j, k), 2.0) + + pow(B_final[2](i, j, k), 2.0); + } + } + } + + return sqrt(energy_final / energy_start); +} + +double FDTD::calc_max_value(Field& field) +{ + double max_val = 0.0; + +#pragma omp parallel for collapse(2) + for (int i = 0; i < parameters.Ni - pml_size_i; i++) + { + for (int j = 0; j < parameters.Nj - pml_size_j; j++) + { + for (int k = 0; k < parameters.Nk - pml_size_k; k++) + { + if (fabs(field(i, j, k)) > fabs(max_val)) + { + max_val = field(i, j, k); + } + } + } + } + return max_val; +} + std::vector> FDTD::update_fields(const int time) { if (time < 0) @@ -310,211 +393,168 @@ std::vector> FDTD::update_fields(const int time) int size_k_main[2] = { 0, parameters.Nk }; for (int t = 0; t < time; t++) { - std::cout << "Iteration: " << t + 1 << std::endl; - - update_B(size_i_main, size_j_main, size_k_main, t); + update_B(size_i_main, size_j_main, size_k_main); update_E(size_i_main, size_j_main, size_k_main, t); - update_B(size_i_main, size_j_main, size_k_main, t); + update_B(size_i_main, size_j_main, size_k_main); std::vector new_iteration{ Ex, Ey, Ez, Bx, By, Bz }; return_data.push_back(new_iteration); } return return_data; } + + // Defining areas of computation + int size_i_main[] = { pml_size_i, parameters.Ni - pml_size_i }; + int size_j_main[] = { pml_size_j, parameters.Nj - pml_size_j }; + int size_k_main[] = { pml_size_k, parameters.Nk - pml_size_k }; - int size_i_main[2] = { pml_size_i, parameters.Ni - pml_size_i }; - int size_j_main[2] = { pml_size_j, parameters.Nj - pml_size_j }; - int size_k_main[2] = { pml_size_k, parameters.Nk - pml_size_k }; + int size_i_solid[] = { 0, parameters.Ni }; + int size_j_solid[] = { 0, parameters.Nj }; + int size_k_solid[] = { 0, parameters.Nk }; - /*int size_i_main[2] = { 0, parameters.Ni }; - int size_j_main[2] = { 0, parameters.Nj }; - int size_k_main[2] = { 0, parameters.Nk };*/ + int size_i_part_from_start[] = { 0, parameters.Ni - pml_size_i }; + int size_i_part_from_end[] = { pml_size_i, parameters.Ni }; - int size_i_solid[2] = { 0, parameters.Ni }; - int size_j_solid[2] = { 0, parameters.Nj }; - int size_k_solid[2] = { 0, parameters.Nk }; + int size_k_part_from_start[] = { 0, parameters.Nk - pml_size_k }; + int size_k_part_from_end[] = { pml_size_k, parameters.Nk }; - int size_xy_lower_k_pml[2] = { 0, pml_size_k }; - int size_xy_upper_k_pml[2] = { parameters.Nk - pml_size_k, parameters.Nk }; + int size_xy_lower_k_pml[] = { 0, pml_size_k }; + int size_xy_upper_k_pml[] = { parameters.Nk - pml_size_k, parameters.Nk }; - int size_yz_lower_i_pml[2] = { 0, pml_size_i }; - int size_yz_upper_i_pml[2] = { parameters.Ni - pml_size_i, parameters.Ni }; + int size_yz_lower_i_pml[] = { 0, pml_size_i }; + int size_yz_upper_i_pml[] = { parameters.Ni - pml_size_i, parameters.Ni }; - int size_zx_lower_j_pml[2] = { 0, pml_size_j }; - int size_zx_upper_j_pml[2] = { parameters.Nj - pml_size_j, parameters.Nj }; + int size_zx_lower_j_pml[] = { 0, pml_size_j }; + int size_zx_upper_j_pml[] = { parameters.Nj - pml_size_j, parameters.Nj }; - std::function calc_distant_i_up = [=](int i, int j, int k) { + // Definition of functions for calculating the distance to the interface + std::function calc_distant_i_up = + [=](int i, int j, int k) { return i + 1 + pml_size_i - parameters.Ni; }; - std::function calc_distant_j_up = [=](int i, int j, int k) { + std::function calc_distant_j_up = + [=](int i, int j, int k) { return j + 1 + pml_size_j - parameters.Nj; }; - std::function calc_distant_k_up = [=](int i, int j, int k) { + std::function calc_distant_k_up = + [=](int i, int j, int k) { return k + 1 + pml_size_k - parameters.Nk; }; - std::function calc_distant_i_low = [=](int i, int j, int k) { + std::function calc_distant_i_low = + [=](int i, int j, int k) { return pml_size_i - i; }; - std::function calc_distant_j_low = [=](int i, int j, int k) { + std::function calc_distant_j_low = + [=](int i, int j, int k) { return pml_size_j - j; }; - std::function calc_distant_k_low = [=](int i, int j, int k) { + std::function calc_distant_k_low = + [=](int i, int j, int k) { return pml_size_k - k; }; - double SGm_Ex = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) * FDTDconst::C * FDTDconst::EPS0 + // Calculation of maximum permittivity + double SGm_x = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) / (static_cast(pml_size_i) * parameters.dx); - double SGm_Ey = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) * FDTDconst::C * FDTDconst::EPS0 + double SGm_y = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) / (static_cast(pml_size_j) * parameters.dy); - double SGm_Ez = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) * FDTDconst::C * FDTDconst::EPS0 + double SGm_z = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) / (static_cast(pml_size_k) * parameters.dz); - double SGm_Bx = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) * FDTDconst::C * FDTDconst::MU0 - / (static_cast(pml_size_i) * parameters.dx); - double SGm_By = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) * FDTDconst::C * FDTDconst::MU0 - / (static_cast(pml_size_j) * parameters.dy); - double SGm_Bz = -(FDTDconst::N + 1.0) / 2.0 * log(FDTDconst::R) * FDTDconst::C * FDTDconst::MU0 - / (static_cast(pml_size_k) * parameters.dz); - - double SG_x[] = { SGm_Ex, SGm_Bx }; - double SG_y[] = { SGm_Ey, SGm_By }; - double SG_z[] = { SGm_Ez, SGm_Bz }; - - set_sigma_z(size_i_solid, size_j_solid, size_xy_lower_k_pml, SG_z, calc_distant_k_low); - set_sigma_y(size_i_solid, size_zx_lower_j_pml, size_k_solid, SG_y, calc_distant_j_low); - set_sigma_x(size_yz_lower_i_pml, size_j_solid, size_k_solid, SG_x, calc_distant_i_low); - - set_sigma_z(size_i_solid, size_j_solid, size_xy_upper_k_pml, SG_z, calc_distant_k_up); - set_sigma_y(size_i_solid, size_zx_upper_j_pml, size_k_solid, SG_y, calc_distant_j_up); - set_sigma_x(size_yz_upper_i_pml, size_j_solid, size_k_solid, SG_x, calc_distant_i_up); - - + // Calculation of permittivity in the cells + set_sigma_z(size_i_solid, size_j_solid, size_xy_lower_k_pml, + SGm_z, calc_distant_k_low); + set_sigma_y(size_i_solid, size_zx_lower_j_pml, size_k_solid, + SGm_y, calc_distant_j_low); + set_sigma_x(size_yz_lower_i_pml, size_j_solid, size_k_solid, + SGm_x, calc_distant_i_low); + + set_sigma_z(size_i_solid, size_j_solid, size_xy_upper_k_pml, + SGm_z, calc_distant_k_up); + set_sigma_y(size_i_solid, size_zx_upper_j_pml, size_k_solid, + SGm_y, calc_distant_j_up); + set_sigma_x(size_yz_upper_i_pml, size_j_solid, size_k_solid, + SGm_x, calc_distant_i_up); + +#ifdef __DEBUG_PML__ double max_val_1 = 0.0; double max_val_2 = 0.0; - - int tmax1 = 0; - int tmax2 = 0; + int t_start = 50; + int t_final = 280; +#endif // __DEBUG_PML__ for (int t = 0; t < time; t++) { std::cout << "Iteration: " << t + 1 << std::endl; - update_B(size_i_main, size_j_main, size_k_main, t); + update_B(size_i_main, size_j_main, size_k_main); - update_B_PML(size_i_solid, size_j_solid, size_xy_lower_k_pml); - update_B_PML(size_i_solid, size_zx_lower_j_pml, size_k_solid); - update_B_PML(size_yz_lower_i_pml, size_j_solid, size_k_solid); + update_B_PML(size_i_part_from_start, size_j_solid, size_xy_lower_k_pml); + update_B_PML(size_i_main, size_zx_lower_j_pml, size_k_main); + update_B_PML(size_yz_lower_i_pml, size_j_solid, size_k_part_from_end); - update_B_PML(size_i_solid, size_j_solid, size_xy_upper_k_pml); - update_B_PML(size_i_solid, size_zx_upper_j_pml, size_k_solid); - update_B_PML(size_yz_upper_i_pml, size_j_solid, size_k_solid); + update_B_PML(size_i_part_from_end, size_j_solid, size_xy_upper_k_pml); + update_B_PML(size_i_main, size_zx_upper_j_pml, size_k_main); + update_B_PML(size_yz_upper_i_pml, size_j_solid, size_k_part_from_start); update_E(size_i_main, size_j_main, size_k_main, t); - update_E_PML(size_i_solid, size_j_solid, size_xy_lower_k_pml); - update_E_PML(size_i_solid, size_zx_lower_j_pml, size_k_solid); - update_E_PML(size_yz_lower_i_pml, size_j_solid, size_k_solid); + update_E_PML(size_i_part_from_start, size_j_solid, size_xy_lower_k_pml); + update_E_PML(size_i_main, size_zx_lower_j_pml, size_k_main); + update_E_PML(size_yz_lower_i_pml, size_j_solid, size_k_part_from_end); - update_E_PML(size_i_solid, size_j_solid, size_xy_upper_k_pml); - update_E_PML(size_i_solid, size_zx_upper_j_pml, size_k_solid); - update_E_PML(size_yz_upper_i_pml, size_j_solid, size_k_solid); + update_E_PML(size_i_part_from_end, size_j_solid, size_xy_upper_k_pml); + update_E_PML(size_i_main, size_zx_upper_j_pml, size_k_main); + update_E_PML(size_yz_upper_i_pml, size_j_solid, size_k_part_from_start); - update_B(size_i_main, size_j_main, size_k_main, t); + update_B(size_i_main, size_j_main, size_k_main); - if (t < 110 && t > 35) +#ifdef __DEBUG_PML__ + if (t == t_start) { - for (int i = 0; i < parameters.Ni - pml_size_i; i++) - { - for (int j = 0; j < parameters.Nj - pml_size_j; j++) - { - for (int k = 0; k < parameters.Nk - pml_size_k; k++) - { - if (fabs(Ex(i, j, k)) > fabs(max_val_1)) - { - max_val_1 = Ex(i, j, k); - tmax1 = t; - } - if (fabs(Ey(i, j, k)) > fabs(max_val_1)) - { - max_val_1 = Ey(i, j, k); - tmax1 = t; - } - if (fabs(Ez(i, j, k)) > fabs(max_val_1)) - { - max_val_1 = Ez(i, j, k); - tmax1 = t; - } - if (fabs(Bx(i, j, k)) > fabs(max_val_1)) - { - max_val_1 = Bx(i, j, k); - tmax1 = t; - } - if (fabs(By(i, j, k)) > fabs(max_val_1)) - { - max_val_1 = By(i, j, k); - tmax1 = t; - } - if (fabs(Bz(i, j, k)) > fabs(max_val_1)) - { - max_val_1 = Bz(i, j, k); - tmax1 = t; - } - } - } - } + max_val_1 = calc_max_value(Ex); } - if (t > 200) + if (t == t_final) { - for (int i = 0; i < parameters.Ni - pml_size_i; i++) - { - for (int j = 0; j < parameters.Nj - pml_size_j; j++) - { - for (int k = 0; k < parameters.Nk - pml_size_k; k++) - { - if (fabs(Ex(i, j, k)) > fabs(max_val_2)) - { - max_val_2 = Ex(i, j, k); - tmax2 = t; - } - if (fabs(Ey(i, j, k)) > fabs(max_val_2)) - { - max_val_2 = Ey(i, j, k); - tmax2 = t; - } - if (fabs(Ez(i, j, k)) > fabs(max_val_2)) - { - max_val_2 = Ez(i, j, k); - tmax2 = t; - } - if (fabs(Bx(i, j, k)) > fabs(max_val_2)) - { - max_val_2 = Bx(i, j, k); - tmax2 = t; - } - if (fabs(By(i, j, k)) > fabs(max_val_2)) - { - max_val_2 = By(i, j, k); - tmax2 = t; - } - if (fabs(Bz(i, j, k)) > fabs(max_val_2)) - { - max_val_2 = Bz(i, j, k); - tmax2 = t; - } - - } - } - } + max_val_2 = calc_max_value(Ex); } - +#ifdef __OUTPUT_SIGMA__ + if (t == 0) break; + std::vector new_iteration{ EsigmaX, EsigmaY, EsigmaZ, + BsigmaX, BsigmaY, BsigmaZ }; +#endif // __OUTPUT_SIGMA__ +#endif //__DEBUG_PML__ + +#ifndef __OUTPUT_SIGMA__ std::vector new_iteration{ Ex, Ey, Ez, Bx, By, Bz }; +#endif // __OUTPUT_SIGMA__ + return_data.push_back(new_iteration); - //if (t == 0) break; } - std::cout << "before: " << max_val_1 << " | " << tmax1 << std::endl; - std::cout << "after: " << max_val_2 << " | " << tmax2 << std::endl; +#ifdef __DEBUG_PML__ + Field E_start[] = { return_data[t_start][0], + return_data[t_start][1], + return_data[t_start][2] }; + Field B_start[] = { return_data[t_start][3], + return_data[t_start][4], + return_data[t_start][5] }; + + Field E_final[] = { return_data[t_final][0], + return_data[t_final][1], + return_data[t_final][2] }; + Field B_final[] = { return_data[t_final][3], + return_data[t_final][4], + return_data[t_final][5] }; + + std::cout << "reflection: " << + calc_reflection(E_start, B_start, E_final, B_final) + << std::endl; + + std::cout << "before: " << max_val_1 << std::endl; + std::cout << "after: " << max_val_2 << std::endl; +#endif // __DEBUG_PML__ return return_data; }