diff --git a/PlotScript/Animations/animation_Bx.gif b/PlotScript/Animations/animation_Bx.gif new file mode 100644 index 0000000..ebb10f6 Binary files /dev/null and b/PlotScript/Animations/animation_Bx.gif differ diff --git a/PlotScript/Animations/animation_Ex.gif b/PlotScript/Animations/animation_Ex.gif index 6a36da6..03d3fcd 100644 Binary files a/PlotScript/Animations/animation_Ex.gif and b/PlotScript/Animations/animation_Ex.gif differ diff --git a/PlotScript/src/Release/sample.exe b/PlotScript/src/Release/sample.exe index 935523b..2e15bb8 100644 Binary files a/PlotScript/src/Release/sample.exe and b/PlotScript/src/Release/sample.exe differ diff --git a/PlotScript/visualization.py b/PlotScript/visualization.py index caa4f0f..853c4d3 100644 --- a/PlotScript/visualization.py +++ b/PlotScript/visualization.py @@ -3,7 +3,6 @@ import pandas as pd import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation -import numpy as np components_num = {"Ex": '1', "Ey": '2', "Ez": '3', "Bx": '4', "By": '5', "Bz": '6'} @@ -19,7 +18,7 @@ def update(frame): ax.clear() data = pd.read_csv(file_names[frame], sep=';') data = data.apply(lambda x: x.str.replace(',', '.').astype(float) if x.dtype == 'object' else x) - ax.imshow(data, cmap='RdBu', vmin=-0.5, vmax=0.5, interpolation='none') + ax.imshow(data, cmap='viridis', vmin=-0.07, vmax=0.07, interpolation='none') # RdBu ax.set_title('Frame {}'.format(frame + 1)) ani = FuncAnimation(fig, update, frames=len(file_names), interval=50) @@ -50,6 +49,5 @@ def execute_cpp(grid_size, iters_num, single_iteration_flag=True): if __name__ == '__main__': - execute_cpp(75, 410, False) + #execute_cpp(75, 410, False) get_animation("Ex") - # get_heatmap("Ex", 86) diff --git a/samples/sample.cpp b/samples/sample.cpp index 7add225..97c5a93 100644 --- a/samples/sample.cpp +++ b/samples/sample.cpp @@ -140,7 +140,9 @@ void spherical_wave(int n, int it, char* base_path = "") // Initialization of the structures and method double d = FDTDconst::C; + double boundary = static_cast(n) / 2.0 * d; + Parameters params { n, @@ -182,7 +184,7 @@ int main(int argc, char* argv[]) { #ifdef __USE_SPHERICAL_WAVE__ int N = 70; - int Iterations = 10; + int Iterations = 410; spherical_wave(N, Iterations, "../../PlotScript/"); #endif diff --git a/src/FDTD.cpp b/src/FDTD.cpp index 45c3b9e..f763286 100644 --- a/src/FDTD.cpp +++ b/src/FDTD.cpp @@ -58,7 +58,7 @@ void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2], 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) / + pow((static_cast(dist(i, j, k))) / static_cast(pml_size_i), FDTDconst::N); } } @@ -78,7 +78,7 @@ void FDTD::set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2], 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) / + pow((static_cast(dist(i, j, k))) / static_cast(pml_size_j), FDTDconst::N); } } @@ -95,11 +95,11 @@ 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 * - pow((static_cast(dist(i, j, k))) - / static_cast(pml_size_k), FDTDconst::N); + pow((static_cast(dist(i, j, k))) / + static_cast(pml_size_k), FDTDconst::N); BsigmaZ(i, j, k) = SGm * - pow((static_cast(dist(i, j, k)) + 0.5) - / static_cast(pml_size_k), FDTDconst::N); + pow((static_cast(dist(i, j, k))) / + static_cast(pml_size_k), FDTDconst::N); } } } @@ -108,12 +108,6 @@ void FDTD::set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2], FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) : parameters(_parameters), dt(_dt), pml_percent(_pml_percent) { - if (parameters.ax > parameters.bx || - parameters.ay > parameters.by || - parameters.az > parameters.bz) - { - throw std::exception("ERROR: invalid parameters"); - } if (parameters.Ni <= 0 || parameters.Nj <= 0 || parameters.Nk <= 0 || @@ -234,6 +228,8 @@ void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) double dy = parameters.dy; double dz = parameters.dz; + double PMLcoef2_x, PMLcoef2_y, PMLcoef2_z; + #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) { @@ -241,33 +237,36 @@ void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) { for (int k = bounds_k[0]; k < bounds_k[1]; k++) { - if (EsigmaX(i, j, k) != 0) - { - 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)) + - (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)) + - (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)) - - (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)) - - (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)) + - (1.0 - PMLcoef(EsigmaZ(i, j, k))) / (EsigmaZ(i, j, k) * dz) * - (Bx(i, j, k) - Bx(i, j, k - 1)); - } + if (EsigmaX(i, j, k) != 0.0) + PMLcoef2_x = (1.0 - PMLcoef(EsigmaX(i, j, k))) / (EsigmaX(i, j, k) * dx); + else + PMLcoef2_x = FDTDconst::C * dt / dx; + + if (EsigmaY(i, j, k) != 0.0) + PMLcoef2_y = (1.0 - PMLcoef(EsigmaY(i, j, k))) / (EsigmaY(i, j, k) * dy); + else + PMLcoef2_y = FDTDconst::C * dt / dy; + + if (EsigmaZ(i, j, k) != 0.0) + PMLcoef2_z = (1.0 - PMLcoef(EsigmaZ(i, j, k))) / (EsigmaZ(i, j, k) * dz); + else + PMLcoef2_z = FDTDconst::C * dt / dz; + + Eyx(i, j, k) = Eyx(i, j, k) * PMLcoef(EsigmaX(i, j, k)) - + PMLcoef2_x * (Bz(i, j, k) - Bz(i - 1, j, k)); + Ezx(i, j, k) = Ezx(i, j, k) * PMLcoef(EsigmaX(i, j, k)) + + PMLcoef2_x * (By(i, j, k) - By(i - 1, j, k)); + + Exy(i, j, k) = Exy(i, j, k) * PMLcoef(EsigmaY(i, j, k)) + + PMLcoef2_y * (Bz(i, j, k) - Bz(i, j - 1, k)); + Ezy(i, j, k) = Ezy(i, j, k) * PMLcoef(EsigmaY(i, j, k)) - + PMLcoef2_y * (Bx(i, j, k) - Bx(i, j - 1, k)); + + Exz(i, j, k) = Exz(i, j, k) * PMLcoef(EsigmaZ(i, j, k)) - + PMLcoef2_z * (By(i, j, k) - By(i, j, k - 1)); + Eyz(i, j, k) = Eyz(i, j, k) * PMLcoef(EsigmaZ(i, j, k)) + + PMLcoef2_z * (Bx(i, j, k) - Bx(i, j, k - 1)); + Ex(i, j, k) = Exz(i, j, k) + Exy(i, j, k); Ey(i, j, k) = Eyx(i, j, k) + Eyz(i, j, k); Ez(i, j, k) = Ezy(i, j, k) + Ezx(i, j, k); @@ -281,6 +280,8 @@ void FDTD::update_B_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; + + double PMLcoef2_x, PMLcoef2_y, PMLcoef2_z; #pragma omp parallel for collapse(2) for (int i = bounds_i[0]; i < bounds_i[1]; i++) @@ -290,32 +291,35 @@ void FDTD::update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2]) for (int k = bounds_k[0]; k < bounds_k[1]; k++) { if (BsigmaX(i, j, k) != 0.0) - { - 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)) - - (1.0 - PMLcoef(BsigmaX(i, j, k))) / (BsigmaX(i, j, k) * dx) * - (Ey(i + 1, j, k) - Ey(i, j, k)); - } + PMLcoef2_x = (1.0 - PMLcoef(BsigmaX(i, j, k))) / (BsigmaX(i, j, k) * dx); + else + PMLcoef2_x = FDTDconst::C * dt / dx; + if (BsigmaY(i, j, k) != 0.0) - { - 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)) + - (1.0 - PMLcoef(BsigmaY(i, j, k))) / (BsigmaY(i, j, k) * dy) * - (Ex(i, j + 1, k) - Ex(i, j, k)); - } + PMLcoef2_y = (1.0 - PMLcoef(BsigmaY(i, j, k))) / (BsigmaY(i, j, k) * dy); + else + PMLcoef2_y = FDTDconst::C * dt / dy; + if (BsigmaZ(i, j, k) != 0.0) - { - 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)) - - (1.0 - PMLcoef(BsigmaZ(i, j, k))) / (BsigmaZ(i, j, k) * dz) * - (Ex(i, j, k + 1) - Ex(i, j, k)); - } + PMLcoef2_z = (1.0 - PMLcoef(BsigmaZ(i, j, k))) / (BsigmaZ(i, j, k) * dz); + else + PMLcoef2_z = FDTDconst::C * dt / dz; + + Byx(i, j, k) = Byx(i, j, k) * PMLcoef(BsigmaX(i, j, k)) + + PMLcoef2_x * (Ez(i + 1, j, k) - Ez(i, j, k)); + Bzx(i, j, k) = Bzx(i, j, k) * PMLcoef(BsigmaX(i, j, k)) - + PMLcoef2_x * (Ey(i + 1, j, k) - Ey(i, j, k)); + + Bxy(i, j, k) = Bxy(i, j, k) * PMLcoef(BsigmaY(i, j, k)) - + PMLcoef2_y * (Ez(i, j + 1, k) - Ez(i, j, k)); + Bzy(i, j, k) = Bzy(i, j, k) * PMLcoef(BsigmaY(i, j, k)) + + PMLcoef2_y * (Ex(i, j + 1, k) - Ex(i, j, k)); + + Bxz(i, j, k) = Bxz(i, j, k) * PMLcoef(BsigmaZ(i, j, k)) + + PMLcoef2_z * (Ey(i, j, k + 1) - Ey(i, j, k)); + Byz(i, j, k) = Byz(i, j, k) * PMLcoef(BsigmaZ(i, j, k)) - + PMLcoef2_z * (Ex(i, j, k + 1) - Ex(i, j, k)); + Bx(i, j, k) = Bxy(i, j, k) + Bxz(i, j, k); By(i, j, k) = Byz(i, j, k) + Byx(i, j, k); Bz(i, j, k) = Bzx(i, j, k) + Bzy(i, j, k); @@ -519,17 +523,23 @@ std::vector> FDTD::update_fields(const int time) { 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__ +#ifdef __OUTPUT_SIGMA__ + + EsigmaX(parameters.Ni / 2, 1, 1) = 100; + + std::vector new_iteration{ EsigmaX, EsigmaY, EsigmaZ, + BsigmaX, BsigmaY, BsigmaZ }; + return_data.push_back(new_iteration); + if (t == 0) break; +#endif // __OUTPUT_SIGMA__ + return_data.push_back(new_iteration); }