pml first implementation

Amazingkivas · Mar 11, 2024 · 5903ef6 · 5903ef6
1 parent 0a60d86
commit 5903ef6
Show file tree

Hide file tree

Showing 88 changed files with 735 additions and 6,883 deletions.
diff --git a/PlotScript/Animations/animation_Bx.gif b/PlotScript/Animations/animation_Bx.gif
diff --git a/PlotScript/Animations/animation_By.gif b/PlotScript/Animations/animation_By.gif
diff --git a/PlotScript/Animations/animation_Bz.gif b/PlotScript/Animations/animation_Bz.gif
diff --git a/PlotScript/Animations/animation_Ex.gif b/PlotScript/Animations/animation_Ex.gif
diff --git a/PlotScript/Animations/animation_Ey.gif b/PlotScript/Animations/animation_Ey.gif
diff --git a/PlotScript/Animations/animation_Ez.gif b/PlotScript/Animations/animation_Ez.gif
diff --git a/PlotScript/OutFiles_1/10.csv b/PlotScript/OutFiles_1/10.csv
diff --git a/PlotScript/OutFiles_1/11.csv b/PlotScript/OutFiles_1/11.csv
diff --git a/PlotScript/OutFiles_1/12.csv b/PlotScript/OutFiles_1/12.csv
diff --git a/PlotScript/OutFiles_1/2.csv b/PlotScript/OutFiles_1/2.csv
diff --git a/PlotScript/OutFiles_1/200.csv b/PlotScript/OutFiles_1/200.csv
diff --git a/PlotScript/OutFiles_1/3.csv b/PlotScript/OutFiles_1/3.csv
diff --git a/PlotScript/OutFiles_1/4.csv b/PlotScript/OutFiles_1/4.csv
diff --git a/PlotScript/OutFiles_1/5.csv b/PlotScript/OutFiles_1/5.csv
diff --git a/PlotScript/OutFiles_1/6.csv b/PlotScript/OutFiles_1/6.csv
diff --git a/PlotScript/OutFiles_1/7.csv b/PlotScript/OutFiles_1/7.csv
diff --git a/PlotScript/OutFiles_1/8.csv b/PlotScript/OutFiles_1/8.csv
diff --git a/PlotScript/OutFiles_1/9.csv b/PlotScript/OutFiles_1/9.csv
diff --git a/PlotScript/OutFiles_2/10.csv b/PlotScript/OutFiles_2/10.csv
diff --git a/PlotScript/OutFiles_2/11.csv b/PlotScript/OutFiles_2/11.csv
diff --git a/PlotScript/OutFiles_2/12.csv b/PlotScript/OutFiles_2/12.csv
diff --git a/PlotScript/OutFiles_2/2.csv b/PlotScript/OutFiles_2/2.csv
diff --git a/PlotScript/OutFiles_2/200.csv b/PlotScript/OutFiles_2/200.csv
diff --git a/PlotScript/OutFiles_2/3.csv b/PlotScript/OutFiles_2/3.csv
diff --git a/PlotScript/OutFiles_2/4.csv b/PlotScript/OutFiles_2/4.csv
diff --git a/PlotScript/OutFiles_2/5.csv b/PlotScript/OutFiles_2/5.csv
diff --git a/PlotScript/OutFiles_2/6.csv b/PlotScript/OutFiles_2/6.csv
diff --git a/PlotScript/OutFiles_2/7.csv b/PlotScript/OutFiles_2/7.csv
diff --git a/PlotScript/OutFiles_2/8.csv b/PlotScript/OutFiles_2/8.csv
diff --git a/PlotScript/OutFiles_2/9.csv b/PlotScript/OutFiles_2/9.csv
diff --git a/PlotScript/OutFiles_3/10.csv b/PlotScript/OutFiles_3/10.csv
diff --git a/PlotScript/OutFiles_3/11.csv b/PlotScript/OutFiles_3/11.csv
diff --git a/PlotScript/OutFiles_3/12.csv b/PlotScript/OutFiles_3/12.csv
diff --git a/PlotScript/OutFiles_3/2.csv b/PlotScript/OutFiles_3/2.csv
diff --git a/PlotScript/OutFiles_3/200.csv b/PlotScript/OutFiles_3/200.csv
diff --git a/PlotScript/OutFiles_3/3.csv b/PlotScript/OutFiles_3/3.csv
diff --git a/PlotScript/OutFiles_3/4.csv b/PlotScript/OutFiles_3/4.csv
diff --git a/PlotScript/OutFiles_3/5.csv b/PlotScript/OutFiles_3/5.csv
diff --git a/PlotScript/OutFiles_3/6.csv b/PlotScript/OutFiles_3/6.csv
diff --git a/PlotScript/OutFiles_3/7.csv b/PlotScript/OutFiles_3/7.csv
diff --git a/PlotScript/OutFiles_3/8.csv b/PlotScript/OutFiles_3/8.csv
diff --git a/PlotScript/OutFiles_3/9.csv b/PlotScript/OutFiles_3/9.csv
diff --git a/PlotScript/OutFiles_4/10.csv b/PlotScript/OutFiles_4/10.csv
diff --git a/PlotScript/OutFiles_4/11.csv b/PlotScript/OutFiles_4/11.csv
diff --git a/PlotScript/OutFiles_4/12.csv b/PlotScript/OutFiles_4/12.csv
diff --git a/PlotScript/OutFiles_4/2.csv b/PlotScript/OutFiles_4/2.csv
diff --git a/PlotScript/OutFiles_4/200.csv b/PlotScript/OutFiles_4/200.csv
diff --git a/PlotScript/OutFiles_4/3.csv b/PlotScript/OutFiles_4/3.csv
diff --git a/PlotScript/OutFiles_4/4.csv b/PlotScript/OutFiles_4/4.csv
diff --git a/PlotScript/OutFiles_4/5.csv b/PlotScript/OutFiles_4/5.csv
diff --git a/PlotScript/OutFiles_4/6.csv b/PlotScript/OutFiles_4/6.csv
diff --git a/PlotScript/OutFiles_4/7.csv b/PlotScript/OutFiles_4/7.csv
diff --git a/PlotScript/OutFiles_4/8.csv b/PlotScript/OutFiles_4/8.csv
diff --git a/PlotScript/OutFiles_4/9.csv b/PlotScript/OutFiles_4/9.csv
diff --git a/PlotScript/OutFiles_5/10.csv b/PlotScript/OutFiles_5/10.csv
diff --git a/PlotScript/OutFiles_5/11.csv b/PlotScript/OutFiles_5/11.csv
diff --git a/PlotScript/OutFiles_5/12.csv b/PlotScript/OutFiles_5/12.csv
diff --git a/PlotScript/OutFiles_5/2.csv b/PlotScript/OutFiles_5/2.csv
diff --git a/PlotScript/OutFiles_5/200.csv b/PlotScript/OutFiles_5/200.csv
diff --git a/PlotScript/OutFiles_5/3.csv b/PlotScript/OutFiles_5/3.csv
diff --git a/PlotScript/OutFiles_5/4.csv b/PlotScript/OutFiles_5/4.csv
diff --git a/PlotScript/OutFiles_5/5.csv b/PlotScript/OutFiles_5/5.csv
diff --git a/PlotScript/OutFiles_5/6.csv b/PlotScript/OutFiles_5/6.csv
diff --git a/PlotScript/OutFiles_5/7.csv b/PlotScript/OutFiles_5/7.csv
diff --git a/PlotScript/OutFiles_5/8.csv b/PlotScript/OutFiles_5/8.csv
diff --git a/PlotScript/OutFiles_5/9.csv b/PlotScript/OutFiles_5/9.csv
diff --git a/PlotScript/OutFiles_6/10.csv b/PlotScript/OutFiles_6/10.csv
diff --git a/PlotScript/OutFiles_6/11.csv b/PlotScript/OutFiles_6/11.csv
diff --git a/PlotScript/OutFiles_6/12.csv b/PlotScript/OutFiles_6/12.csv
diff --git a/PlotScript/OutFiles_6/2.csv b/PlotScript/OutFiles_6/2.csv
diff --git a/PlotScript/OutFiles_6/200.csv b/PlotScript/OutFiles_6/200.csv
diff --git a/PlotScript/OutFiles_6/3.csv b/PlotScript/OutFiles_6/3.csv
diff --git a/PlotScript/OutFiles_6/4.csv b/PlotScript/OutFiles_6/4.csv
diff --git a/PlotScript/OutFiles_6/5.csv b/PlotScript/OutFiles_6/5.csv
diff --git a/PlotScript/OutFiles_6/6.csv b/PlotScript/OutFiles_6/6.csv
diff --git a/PlotScript/OutFiles_6/7.csv b/PlotScript/OutFiles_6/7.csv
diff --git a/PlotScript/OutFiles_6/8.csv b/PlotScript/OutFiles_6/8.csv
diff --git a/PlotScript/OutFiles_6/9.csv b/PlotScript/OutFiles_6/9.csv
diff --git a/PlotScript/src/Release/sample.exe b/PlotScript/src/Release/sample.exe
diff --git a/PlotScript/visualization.py b/PlotScript/visualization.py
@@ -12,8 +12,6 @@ def get_animation(component):
     file_names_base = sorted(glob.glob(folder_path + '/*.csv'))
     file_names = sorted(file_names_base, key=lambda x: int(x.split('\\')[-1].split('.')[0]))
 
-    data_frames = [pd.read_csv(file) for file in file_names]
-
     fig, ax = plt.subplots()
 
     def update(frame):
@@ -23,14 +21,25 @@ def update(frame):
         ax.imshow(data, cmap='viridis', interpolation='nearest')
         ax.set_title('Frame {}'.format(frame + 1))
 
-    ani = FuncAnimation(fig, update, frames=len(file_names), interval=20)
+    ani = FuncAnimation(fig, update, frames=len(file_names), interval=50)
+
     plt.show()
     ani.save(f'animations/animation_{component}.gif', writer='imagemagick')
 
 
-def execute_cpp(grid_size, iters_num):
+def get_heatmap(component, iteration):
+    data = pd.read_csv(f'OutFiles_{components_num[component]}/{iteration}.csv', sep=';')
+
+    data = data.applymap(lambda x: float(x.replace(',', '.')) if isinstance(x, str) else x)
+
+    plt.imshow(data, cmap='viridis', interpolation='nearest')
+    plt.colorbar()
+    plt.show()
+
+
+def execute_cpp(grid_size, iters_num, single_iteration_flag=True):
     cpp_executable = 'src/Release/sample.exe'
-    args = [cpp_executable, str(grid_size), str(iters_num)]
+    args = [cpp_executable, str(grid_size), str(iters_num), str(int(single_iteration_flag))]
     try:
         subprocess.run(args, check=True)
     except subprocess.CalledProcessError:
@@ -40,5 +49,6 @@ def execute_cpp(grid_size, iters_num):
 
 
 if __name__ == '__main__':
-    execute_cpp(100, 12)
-    get_animation("By")
+    #execute_cpp(50, 100, False)
+    get_animation("Ey")
+    #get_heatmap("Bx", 100)
diff --git a/include/FDTD.h b/include/FDTD.h
@@ -4,6 +4,7 @@
 
 #include <vector>
 #include <omp.h>
+#include <functional>
 
 #include "Structures.h"
 
@@ -30,17 +31,39 @@ class FDTD
 {
 private:
     Field Ex, Ey, Ez, Bx, By, Bz;
+
+    Field Exy, Exz, Eyx, Eyz, Ezx, Ezy;
+    Field Bxy, Bxz, Byx, Byz, Bzx, Bzy;
+
     std::vector<Field> Jx;
     std::vector<Field> Jy;
     std::vector<Field> Jz;
+
+    std::vector<double> EsigmaX;
+    std::vector<double> EsigmaY;
+    std::vector<double> EsigmaZ;
+    std::vector<double> BsigmaX;
+    std::vector<double> BsigmaY;
+    std::vector<double> BsigmaZ;
+
     Parameters parameters;
+    PMLparameters pml;
+    double pml_percent;
+    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);
+
+    double PMLcoef(std::vector<double> conductivity, double constant, int distance);
+    void update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t, std::function<int(int, int, int)> dist);
+    void update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t, std::function<int(int, int, int)> dist);
+
 public:
-    FDTD(Parameters _parameters, double _dt);
+    FDTD(Parameters _parameters, double _dt, double _pml_percent);
 
     Field& get_field(Component);
     std::vector<Field>& get_current(Component);
 
-    void update_fields(const int time);
+    std::vector<std::vector<Field>> update_fields(const int time);
 };
diff --git a/include/Structures.h b/include/Structures.h
@@ -4,7 +4,12 @@
 
 namespace FDTDconst
 {
-	const double C = 1.0;  // light speed
+	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 N = 4.0;
 }
 
 namespace FDTDstruct
@@ -22,9 +27,19 @@ namespace FDTDstruct
 		int period;
 		int m;
 		double dt;
-		int iterations = static_cast<int>(static_cast<double>(period) / dt);;
+		int iterations;
 		double period_x = static_cast<double>(m) * FDTDconst::C;
 		double period_y = static_cast<double>(m) * FDTDconst::C;
+		double period_z = static_cast<double>(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 {

diff --git a/include/Writer.h b/include/Writer.h
@@ -72,29 +72,80 @@ void write_z(Field& this_field, std::ofstream& fout)
     fout << std::endl << std::endl;
 }
 
-void write_spherical(FDTD& test, Axis axis, int iteration, std::string base_path)
+void write_spherical(std::vector<std::vector<Field>> data, Axis axis, int iteration, std::string base_path)
 {
     std::ofstream test_fout;
 
-    for (int i = static_cast<int>(Component::EX); i <= static_cast<int>(Component::BZ); ++i)
+    int it = 0;
+    for (std::vector<Field> fields : data)
     {
-        if (axis == Axis::X)
-        {
-            test_fout.open(base_path + "OutFiles_" + std::to_string(i + 1) + "/" + std::to_string(iteration) + ".csv");
-            write_x(test.get_field(static_cast<Component>(i)), test_fout);
-            test_fout.close();
-        }
-        else if (axis == Axis::Y)
+        it++;
+        for (int c = static_cast<int>(Component::EX); c <= static_cast<int>(Component::BZ); ++c)
         {
-            test_fout.open(base_path + "OutFiles_" + std::to_string(i + 1) + "/" + std::to_string(iteration) + ".csv");
-            write_y(test.get_field(static_cast<Component>(i)), test_fout);
-            test_fout.close();
-        }
-        else
-        {
-            test_fout.open(base_path + "OutFiles_" + std::to_string(i + 1) + "/" + std::to_string(iteration) + ".csv");
-            write_z(test.get_field(static_cast<Component>(i)), test_fout);
-            test_fout.close();
+            std::cout << "Write : " << it << " -- " << c << std::endl;
+            if (axis == Axis::X)
+            {
+                test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv");
+                Field field = fields[c];
+                for (int j = 0; j < field.get_Nj(); ++j)
+                {
+                    for (int k = 0; k < field.get_Nk(); ++k)
+                    {
+                        test_fout << field(field.get_Ni() / 2, j, k);
+                        if (k == field.get_Nk() - 1)
+                        {
+                            test_fout << std::endl;
+                        }
+                        else
+                        {
+                            test_fout << ";";
+                        }
+                    }
+                }
+                test_fout.close();
+            }
+            else if (axis == Axis::Y)
+            {
+                test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv");
+                Field field = fields[c];
+                for (int i = 0; i < field.get_Nj(); ++i)
+                {
+                    for (int k = 0; k < field.get_Nk(); ++k)
+                    {
+                        test_fout << field(i, field.get_Nj() / 2, k);
+                        if (k == field.get_Nk() - 1)
+                        {
+                            test_fout << std::endl;
+                        }
+                        else
+                        {
+                            test_fout << ";";
+                        }
+                    }
+                }
+                test_fout.close();
+            }
+            else
+            {
+                test_fout.open(base_path + "OutFiles_" + std::to_string(c + 1) + "/" + std::to_string(it) + ".csv");
+                Field field = fields[c];
+                for (int i = 0; i < field.get_Nj(); ++i)
+                {
+                    for (int j = 0; j < field.get_Nk(); ++j)
+                    {
+                        test_fout << field(i, j, field.get_Nk() / 2);
+                        if (j == field.get_Nk() - 1)
+                        {
+                            test_fout << std::endl;
+                        }
+                        else
+                        {
+                            test_fout << ";";
+                        }
+                    }
+                }
+                test_fout.close();
+            }
         }
     }
 }

diff --git a/include/test_FDTD.h b/include/test_FDTD.h
@@ -20,11 +20,12 @@ class Test_FDTD
 public:
 	Test_FDTD(Parameters);
 
+	void initiialize_current(FDTD& _test, CurrentParameters, int iters,
+		std::function<double(double, double, double, double)>& init_function);
+
 	void initial_filling(FDTD& _test, SelectedFields, int iters,
 		std::function<double(double, double[2])>& init_function);
 
-	void initiialize_current(FDTD& _test, CurrentParameters, int iters, std::function<double(double, double, double)>& init_function);
-
 	double get_max_abs_error(Field& this_field, Component field,
 		std::function<double(double, double, double[2])>& true_function, double time);
 };
diff --git a/samples/sample.cpp b/samples/sample.cpp
@@ -47,7 +47,6 @@ void plane_wave(std::vector<char*> arguments, std::vector<double> numbers, char*
     Component flds_selected[2] = { selected_field_1, selected_field_2 };
 
     // Saving data from the txt file
-
     int grid_sizes[3] = { numbers[0], numbers[1], numbers[2] };   // { Nx, Ny , Nz }
     double sizes_x[2] = { numbers[3], numbers[4] };               // { ax, bx }
     double sizes_y[2] = { numbers[5], numbers[6] };               // { ay, by }
@@ -90,7 +89,7 @@ void plane_wave(std::vector<char*> arguments, std::vector<double> numbers, char*
         (sizes_y[1] - sizes_y[0]) / static_cast<double>(grid_sizes[1]),
         (sizes_z[1] - sizes_z[0]) / static_cast<double>(grid_sizes[2])
     };
-    FDTD method(params, time_step);
+    FDTD method(params, time_step, 0.0);
 
     // Meaningful calculations
     Test_FDTD test(params);
@@ -132,37 +131,42 @@ void spherical_wave(int n, int it, char* base_path = "")
     double T = cur_param.period;
     double Tx = cur_param.period_x;
     double Ty = cur_param.period_y;
-    std::function<double(double, double, double)> cur_func = [T, Tx, Ty](double x, double y, 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);
-    };
-
+    double Tz = cur_param.period_z;
+    double omega1 = 1.0;
+    double omega2 = 0.5;
+    cur_param.iterations = static_cast<int>(static_cast<double>(cur_param.period) / cur_param.dt); //2.0 * M_PI / (omega2 * cur_param.dt)
+    std::function<double(double, double, double, double)> cur_func = [omega1, omega2, 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);
+    }; // sin(omega1 * t) * pow(sin(omega2 * t), 2.0)
+
     // Initialization of the structures and method
     double d = FDTDconst::C;
+    double boundary = static_cast<double>(n) / 2.0 * d;
     Parameters params
     {
         n,
         n,
-        1,
-        -static_cast<double>(n) / 2.0 * d,
-        static_cast<double>(n) / 2.0 * d,
-        -static_cast<double>(n) / 2.0 * d,
-        static_cast<double>(n) / 2.0 * d,
-        0.0,
-        0.0,
+        n,
+        -boundary,
+        boundary,
+        -boundary,
+        boundary,
+        -boundary,
+        boundary,
         d,
         d,
         d
     };
-    FDTD method(params, cur_param.dt);
+    FDTD method(params, cur_param.dt, 0.1);
 
     // Meaningful calculations
     Test_FDTD test(params);
     test.initiialize_current(method, cur_param, it, cur_func);
-    method.update_fields(it);
+    std::vector<std::vector<Field>> data = method.update_fields(it);
     try
     {
-        write_spherical(method, Axis::X, it, base_path);
+        write_spherical(data, Axis::X, it, base_path);
     }
     catch (const std::exception& except)
     {
@@ -179,13 +183,9 @@ int main(int argc, char* argv[])
     if (argc == 1) 
     {
 #ifdef __USE_SPHERICAL_WAVE__
-        int N = 100;
-        int Iterations = 11;
-        for (int I = 2; I <= Iterations; I++)
-        {
-            std::cout << "Iteration: " << I << std::endl;
-            spherical_wave(N, I, "../../PlotScript/");
-        }
+        int N = 50;
+        int Iterations = 200;
+        spherical_wave(N, Iterations, "../../PlotScript/");
 #endif
 
 #ifndef __USE_SPHERICAL_WAVE__
@@ -223,15 +223,11 @@ int main(int argc, char* argv[])
         plane_wave(arguments, numbers, outfile_path);
 #endif
     }
-    else if (argc == 3)
+    else if (argc == 4)
     {
         int N = std::atoi(arguments[1]);
         int Iterations = std::atoi(arguments[2]);
-        for (int I = 2; I <= Iterations; I++)
-        {
-            std::cout << "Iteration: " << I << std::endl;
-            spherical_wave(N, I);
-        }
+        spherical_wave(N, Iterations);
     }
     else
     {