fix memory usage with current

Amazingkivas · Oct 11, 2024 · 13e9022 · 13e9022
1 parent 181f783
commit 13e9022
Show file tree

Hide file tree

Showing 8 changed files with 82 additions and 42 deletions.
diff --git a/PlotScript/Animations/animation_Bx.gif b/PlotScript/Animations/animation_Bx.gif
diff --git a/include/Current.h b/include/Current.h
@@ -0,0 +1,20 @@
+#pragma once
+
+#include <vector>
+
+#include "Field.h"
+
+class Current
+{
+private:
+	std::vector<Field> J;
+	int iterations;
+	int Ni, Nj, Nk;
+	double default_value = 0.0;
+public:
+	Current(int iters = 1, int _Ni = 1, int _Nj = 1, int _Nk = 1);
+
+	Field& operator[] (int iteration);
+
+	double& operator() (int iteration, int i, int j, int k);
+};
diff --git a/include/FDTD.h b/include/FDTD.h
@@ -8,6 +8,7 @@
 
 #include "Writer.h"
 #include "Field.h"
+#include "Current.h"
 
 using namespace FDTDstruct;
 
@@ -16,12 +17,14 @@ class FDTD
 private:
     Field Ex, Ey, Ez, Bx, By, Bz;
 
-    std::vector<Field> Jx;
-    std::vector<Field> Jy;
-    std::vector<Field> Jz;
+    Current Jx;
+    Current Jy;
+    Current Jz;
 
     Field Exy, Exz, Eyx, Eyz, Ezx, Ezy;
     Field Bxy, Bxz, Byx, Byz, Bzx, Bzy;
+
+    // Permittivity and permeability of the medium
     Field EsigmaX, EsigmaY, EsigmaZ, 
           BsigmaX, BsigmaY, BsigmaZ;
 
@@ -44,10 +47,10 @@ class FDTD
         double SGm, std::function<int(int, int, int)> dist);
 
 public:
-    FDTD(Parameters _parameters, double _dt, double _pml_percent);
+    FDTD(Parameters _parameters, double _dt, double _pml_percent, int current_iters = 1);
 
     Field& get_field(Component);
-    std::vector<Field>& get_current(Component);
+    Current& get_current(Component);
 
     std::vector<Field> update_fields(const int time, bool write_result = false, 
         Axis write_axis = Axis::X, std::string base_path = "");

diff --git a/samples/sample.cpp b/samples/sample.cpp
@@ -76,7 +76,7 @@ void spherical_wave(int n, int it, char* base_path = "")
         d,
         d
     };
-    FDTD method(params, cur_param.dt, 0.1);
+    FDTD method(params, cur_param.dt, 0.1, cur_param.iterations);
 
     // Meaningful calculations
     Test_FDTD test(params);

diff --git a/sln/FDTD/CMakeLists.txt b/sln/FDTD/CMakeLists.txt
@@ -1,3 +1,3 @@
 set(target "FDTD")
 
-add_library(${target} STATIC ../../include/FDTD.h ../../include/Field.h ../../src/FDTD.cpp ../../src/Field.cpp)
+add_library(${target} STATIC ../../include/FDTD.h ../../include/Field.h ../../include/Current.h ../../src/FDTD.cpp ../../src/Field.cpp ../../src/Current.cpp)
diff --git a/src/Current.cpp b/src/Current.cpp
@@ -0,0 +1,31 @@
+#include "Current.h"
+
+Current::Current(int iters, int _Ni, int _Nj, int _Nk) 
+	: iterations(iters), Ni(_Ni), Nj(_Nj), Nk(_Nk)
+{
+	J = std::vector<Field>(iterations, Field(Ni, Nj, Nk));
+}
+
+Field& Current::operator[] (int iteration)
+{
+	if (iteration < iterations)
+	{
+		return J[iteration];
+	}
+	else
+	{
+		return Field(Ni, Nj, Nk);
+	}
+}
+
+double& Current::operator() (int iteration, int i, int j, int k)
+{
+	if (iteration < iterations)
+	{
+		return J[iteration](i, j, k);
+	}
+	else
+	{
+		return default_value;
+	}
+}
diff --git a/src/FDTD.cpp b/src/FDTD.cpp
@@ -4,12 +4,10 @@
 #include <exception>
 #include <cmath>
 
-#include "Field.h"
-
 void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2],
     double SGm, std::function<int(int, int, int)> dist)
 {
-#pragma omp parallel for collapse(2)
+#pragma omp parallel for collapse(3)
     for (int i = bounds_i[0]; i < bounds_i[1]; i++)
     {
         for (int j = bounds_j[0]; j < bounds_j[1]; j++)
@@ -29,7 +27,7 @@ void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2],
 void FDTD::set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2],
     double SGm, std::function<int(int, int, int)> dist)
 {
-#pragma omp parallel for collapse(2)
+#pragma omp parallel for collapse(3)
     for (int i = bounds_i[0]; i < bounds_i[1]; i++)
     {
         for (int j = bounds_j[0]; j < bounds_j[1]; j++)
@@ -49,7 +47,7 @@ void FDTD::set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2],
 void FDTD::set_sigma_z(int bounds_i[2], int bounds_j[2], int bounds_k[2],
     double SGm, std::function<int(int, int, int)> dist)
 {
-#pragma omp parallel for collapse(2)
+#pragma omp parallel for collapse(3)
     for (int i = bounds_i[0]; i < bounds_i[1]; i++)
     {
         for (int j = bounds_j[0]; j < bounds_j[1]; j++)
@@ -67,7 +65,7 @@ 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) :
+FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent, int current_iters) :
     parameters(_parameters), dt(_dt), pml_percent(_pml_percent)
 {
     if (parameters.Ni <= 0 ||
@@ -77,6 +75,7 @@ FDTD::FDTD(Parameters _parameters, double _dt, double _pml_percent) :
     {
         throw std::invalid_argument("ERROR: invalid parameters");
     }
+    Jx = Jy = Jz = Current(current_iters, 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
@@ -111,7 +110,7 @@ Field& FDTD::get_field(Component this_field)
     }
 }
 
-std::vector<Field>& FDTD::get_current(Component this_current)
+Current& FDTD::get_current(Component this_current)
 {
     switch (this_current)
     {
@@ -131,20 +130,20 @@ void FDTD::update_E(int bounds_i[2], int bounds_j[2], int bounds_k[2], int t)
     double dy = parameters.dy;
     double dz = parameters.dz;
 
-#pragma omp parallel for collapse(2)
+#pragma omp parallel for collapse(3)
     for (int i = bounds_i[0]; i < bounds_i[1]; i++)
     {
         for (int j = bounds_j[0]; j < bounds_j[1]; j++)
         {
             for (int k = bounds_k[0]; k < bounds_k[1]; k++)
             {
-                Ex(i, j, k) = Ex(i, j, k) - 4.0 * FDTDconst::PI * dt * Jx[t](i, j, k) +
+                Ex(i, j, k) = Ex(i, j, k) - 4.0 * FDTDconst::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 * FDTDconst::PI * dt * Jy[t](i, j, k) +
+                Ey(i, j, k) = Ey(i, j, k) - 4.0 * FDTDconst::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 * FDTDconst::PI * dt * Jz[t](i, j, k) +
+                Ez(i, j, k) = Ez(i, j, k) - 4.0 * FDTDconst::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);
             }
@@ -158,7 +157,7 @@ void FDTD::update_B(int bounds_i[2], int bounds_j[2], int bounds_k[2])
     double dy = parameters.dy;
     double dz = parameters.dz;
 
-#pragma omp parallel for collapse(2)
+#pragma omp parallel for collapse(3)
     for (int i = bounds_i[0]; i < bounds_i[1]; i++)
     {
         for (int j = bounds_j[0]; j < bounds_j[1]; j++)
@@ -192,7 +191,7 @@ void FDTD::update_E_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2])
 
     double PMLcoef2_x, PMLcoef2_y, PMLcoef2_z;
 
-#pragma omp parallel for collapse(2)
+#pragma omp parallel for collapse(3)
     for (int i = bounds_i[0]; i < bounds_i[1]; i++)
     {
         for (int j = bounds_j[0]; j < bounds_j[1]; j++)
@@ -245,7 +244,7 @@ void FDTD::update_B_PML(int bounds_i[2], int bounds_j[2], int bounds_k[2])
 
     double PMLcoef2_x, PMLcoef2_y, PMLcoef2_z;
 
-#pragma omp parallel for collapse(2)
+#pragma omp parallel for collapse(3)
     for (int i = bounds_i[0]; i < bounds_i[1]; i++)
     {
         for (int j = bounds_j[0]; j < bounds_j[1]; j++)
@@ -368,15 +367,15 @@ std::vector<Field> FDTD::update_fields(const int time, bool write_result, Axis w
         return pml_size_k - k;
     };
 
-    // Calculation of maximum permittivity
+    // Calculation of maximum permittivity and permeability
     double SGm_x = -(FDTDconst::N + 1.0) / 2.0 * std::log(FDTDconst::R)
         / (static_cast<double>(pml_size_i) * parameters.dx);
     double SGm_y = -(FDTDconst::N + 1.0) / 2.0 * std::log(FDTDconst::R)
         / (static_cast<double>(pml_size_j) * parameters.dy);
     double SGm_z = -(FDTDconst::N + 1.0) / 2.0 * std::log(FDTDconst::R)
         / (static_cast<double>(pml_size_k) * parameters.dz);
 
-    // Calculation of permittivity in the cells
+    // Calculation of permittivity and permeability 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,

diff --git a/src/test_FDTD.cpp b/src/test_FDTD.cpp
@@ -18,16 +18,10 @@ void Test_FDTD::initiialize_current(FDTD& _test, CurrentParameters cParams, int
 	int max_j = std::floor((Ty / 4.0 - parameters.ay) / parameters.dy);
 	int max_k = std::floor((Tz / 4.0 - parameters.az) / parameters.dz);
 
-	for (int n = 0; n < std::max(iters, cParams.iterations); n++)
-	{
-		_test.get_current(Component::JX).push_back(Field(parameters.Ni, parameters.Nj, parameters.Nk));
-		_test.get_current(Component::JY).push_back(Field(parameters.Ni, parameters.Nj, parameters.Nk));
-		_test.get_current(Component::JZ).push_back(Field(parameters.Ni, parameters.Nj, parameters.Nk));
-	}
-	for (int iter = 1; iter < cParams.iterations; iter++)
+	for (int iter = 0; iter < cParams.iterations; iter++)
 	{
 		Field J(parameters.Ni, parameters.Nj, parameters.Nk);
-
+#pragma omp parallel for collapse(3)
 		for (int i = start_i; i <= max_i; i++)
 		{
 			for (int j = start_j; j <= max_j; j++)
@@ -37,26 +31,19 @@ void Test_FDTD::initiialize_current(FDTD& _test, CurrentParameters cParams, int
 					J(i, j, k) = init_function(static_cast<double>(i) * parameters.dx,
 						static_cast<double>(j) * parameters.dy,
 						static_cast<double>(k) * parameters.dz,
-						static_cast<double>(iter) * cParams.dt);
+						static_cast<double>(iter + 1) * cParams.dt);
 				}
 			}
 		}
-		_test.get_current(Component::JX)[iter - 1] = J;
-		_test.get_current(Component::JY)[iter - 1] = J;
-		_test.get_current(Component::JZ)[iter - 1] = J;
+		_test.get_current(Component::JX)[iter] = J;
+		_test.get_current(Component::JY)[iter] = J;
+		_test.get_current(Component::JZ)[iter] = J;
 	}
 }
 
 void Test_FDTD::initial_filling(FDTD& _test, SelectedFields fields, int iters,
 	std::function<double(double, double[2])>& init_function)
 {
-	for (int n = 0; n < iters; n++)
-	{
-		_test.get_current(Component::JX).push_back(Field(parameters.Ni, parameters.Nj, parameters.Nk));
-		_test.get_current(Component::JY).push_back(Field(parameters.Ni, parameters.Nj, parameters.Nk));
-		_test.get_current(Component::JZ).push_back(Field(parameters.Ni, parameters.Nj, parameters.Nk));
-	}
-
 	set_axis(fields.selected_E, fields.selected_B);
 	set_sign(fields.selected_E, fields.selected_B);
 	if (axis == Axis::X)