Skip to content

Commit

Permalink
fix boundary conditions & convergence study
Browse files Browse the repository at this point in the history
  • Loading branch information
@Amazingkivas
Amazingkivas committed Dec 5, 2023
1 parent 60f96a8 commit 09ecacf
Show file tree
Hide file tree
Showing 7 changed files with 276 additions and 983 deletions.
960 changes: 192 additions & 768 deletions PlotScript/OutFile.csv
View file

Large diffs are not rendered by default.

7 changes: 3 additions & 4 deletions PlotScript/Source.txt
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,7 @@
4.0
4.0
32
32
0.0
1.0
0.0
1.0
2e-13
1e-11
32
9 changes: 2 additions & 7 deletions include/FDTD.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,13 +39,8 @@ class FDTD

Field& get_field(Component);

void periodical_boundary_conditions_E();
void periodical_boundary_conditions_B();
void periodical_boundary_conditions_shiftedE();
void periodical_boundary_conditions_shiftedB();

void update_field(const double&);
void shifted_update_field(const double&);
void update_field(const int&);
void shifted_update_field(const int&);

int get_Ni() { return Ni; }
int get_Nj() { return Nj; }
Expand Down
8 changes: 4 additions & 4 deletions include/test_FDTD.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,19 +17,19 @@ class Test_FDTD
bool shifted;
double max_abs_error = 0.0;

void initial_filling(FDTD& _test, Component field_1, Component field_2, double size_d, double size_wave[2],
void initial_filling(FDTD& _test, Component fields[2], double size_d, double size_wave[2],
std::function<double(double, double[2])>& _init_function);

void start(FDTD& _test, double _t);
void start_test(FDTD& _test, double _t);

double get_shift(Component _field, double step);

void get_max_abs_error(Field& this_field, Component field, double _size_d[2], double size_wave[2],
std::function<double(double, double, double[2])>& _true_function, double _t);

public:
Test_FDTD(FDTD& test, Component field_1, Component field_2, Component field_3,
double size_x[2], double size_y[2], double size_d[2], double t,
Test_FDTD(FDTD& test, Component fields[2], Component field_3,
double size_x[2], double size_y[2], double size_d[2], double time, int iters,
std::function<double(double, double[2])>& init_function,
std::function<double(double, double, double[2])>& true_function, bool _shifted);

Expand Down
51 changes: 31 additions & 20 deletions samples/sample.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,15 +72,6 @@ void write_all(FDTD& test, char axis)
test_fout.close();
}

double func_1(double x, double size[2])
{
return sin(2.0 * M_PI * (x - size[0]) / (size[1] - size[0]));
}
double func_2(double x, double t, double size[2])
{
return sin(2.0 * M_PI * (x - size[0] - FDTD_Const::C * t) / (size[1] - size[0]));
}

#ifndef __DEBUG__
int main(int argc, char* argv[])
#else
Expand All @@ -103,7 +94,7 @@ int main()
}

#ifdef __DEBUG__
const char* argv[5] = {"0", "0", "5", "0", "1"};
const char* argv[5] = {"1", "2", "3", "2", "1"};
#endif

double number;
Expand All @@ -117,32 +108,52 @@ int main()
double arr_y[2] = { numbers[4], numbers[5] };
double arr_d[2] = { (arr_x[1] - arr_x[0]) / arr_N[0], (arr_y[1] - arr_y[0]) / arr_N[1] };

FDTD test_1(arr_N, arr_x, arr_y, numbers[6]);
double courant_condtition = 0.25;
double dt;
double t;

Component fld_1 = static_cast<Component>(std::atoi(argv[1]));
Component fld_2 = static_cast<Component>(std::atoi(argv[2]));
Component fld_3 = static_cast<Component>(std::atoi(argv[3]));
Component flds_selected[2] = { static_cast<Component>(std::atoi(argv[1])), static_cast<Component>(std::atoi(argv[2])) };
Component fld_tested = static_cast<Component>(std::atoi(argv[3]));
bool shifted_flag = static_cast<bool>(std::atoi(argv[4]));

std::function<double(double, double[2])> initial_func = func_1;
std::function<double(double, double, double[2])> true_func = func_2;

Test_FDTD test(test_1, fld_1, fld_2, fld_3, arr_x, arr_y, arr_d, numbers[7], initial_func, true_func, static_cast<bool>(std::atoi(argv[4])));
std::cout << test.get_max_abs_error() << std::endl;
std::function<double(double, double[2])> initial_func =
[](double x, double size[2]) { return sin(2.0 * M_PI * (x - size[0]) / (size[1] - size[0])); };

std::function<double(double, double, double[2])> true_func =
[](double x, double t, double size[2]) { return sin(2.0 * M_PI * (x - size[0] - FDTD_Const::C * t) / (size[1] - size[0])); };


char selected_axis;
if (fld_1 == Component::EY && fld_2 == Component::BZ || fld_1 == Component::EZ && fld_2 == Component::BY)
if (flds_selected[0] == Component::EY && flds_selected[1] == Component::BZ ||
flds_selected[0] == Component::EZ && flds_selected[1] == Component::BY)
{
selected_axis = 'x';
dt = courant_condtition * arr_d[0] / FDTD_Const::C;
t = courant_condtition * (arr_x[1] - arr_x[0]) / FDTD_Const::C;
}
else if (fld_1 == Component::EX && fld_2 == Component::BZ || fld_1 == Component::EZ && fld_2 == Component::BX)
else if (flds_selected[0] == Component::EX && flds_selected[1] == Component::BZ ||
flds_selected[0] == Component::EZ && flds_selected[1] == Component::BX)
{
selected_axis = 'y';
dt = courant_condtition * arr_d[1] / FDTD_Const::C;
t = courant_condtition * (arr_y[1] - arr_y[0]) / FDTD_Const::C;
}
else
{
std::cout << "ERROR" << std::endl;
exit(1);
}


FDTD test_1(arr_N, arr_x, arr_y, dt);


int iter_nums = static_cast<int>(numbers[6]);
Test_FDTD test(test_1, flds_selected, fld_tested, arr_x, arr_y, arr_d, t, iter_nums, initial_func, true_func, shifted_flag);
std::cout << test.get_max_abs_error() << std::endl;


write_all(test_1, selected_axis);

source_fin.close();
Expand Down
180 changes: 23 additions & 157 deletions src/FDTD.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,16 @@ Field& Field::operator= (const Field& other)
return *this;
}

double& Field::operator() (int _i, int _j)
double& Field::operator() (int i, int j)
{
int index = _j + _i * Nj;
int i_isMinusOne = (i == -1);
int j_isMinusOne = (j == -1);
int i_isNi = (i == Ni);
int j_isNj = (j == Nj);
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 index = truly_j + truly_i * Nj;
return field[index];
}

Expand Down Expand Up @@ -55,61 +62,57 @@ Field& FDTD::get_field(Component this_field)
}
}

void FDTD::update_field(const double& time)
void FDTD::update_field(const int& time)
{
for (double t = 0; t < time; t += dt)
for (double t = 0; t <= time; ++t)
{
#pragma omp parallel for
for (int j = 1; j < Nj - 1; ++j)
for (int j = 0; j < Nj; ++j)
{
#pragma ivdep
for (int i = 1; i < Ni - 1; ++i)
for (int i = 0; i < Ni; ++i)
{
Ex(i, j) += FDTD_Const::C * dt * (Bz(i, j + 1) - Bz(i, j - 1)) / (2.0 * dy);
Ey(i, j) -= FDTD_Const::C * dt * (Bz(i + 1, j) - Bz(i - 1, j)) / (2.0 * dx);
Ez(i, j) += FDTD_Const::C * dt * ((By(i + 1, j) - By(i - 1, j)) / (2.0 * dx) - (Bx(i, j + 1) - Bx(i, j - 1)) / (2.0 * dy));
}
}
periodical_boundary_conditions_E();

#pragma omp parallel for
for (int j = 1; j < Nj - 1; ++j)
for (int j = 0; j < Nj; ++j)
{
#pragma ivdep
for (int i = 1; i < Ni - 1; ++i)
for (int i = 0; i < Ni; ++i)
{
Bx(i, j) -= FDTD_Const::C * dt * (Ez(i, j + 1) - Ez(i, j - 1)) / (2.0 * dy);
By(i, j) += FDTD_Const::C * dt * (Ez(i + 1, j) - Ez(i - 1, j)) / (2.0 * dx);
Bz(i, j) -= FDTD_Const::C * dt * ((Ey(i + 1, j) - Ey(i - 1, j)) / (2.0 * dx) - (Ex(i, j + 1) - Ex(i, j - 1)) / (2.0 * dy));
}
}
periodical_boundary_conditions_B();
}
}

void FDTD::shifted_update_field(const double& time)
void FDTD::shifted_update_field(const int& time)
{
for (double t = 0; t <= time; t += dt)
for (double t = 0; t <= time; ++t)
{
#pragma omp parallel for
for (int j = 0; j < Nj - 1; ++j)
for (int j = 0; j < Nj; ++j)
{
#pragma ivdep
for (int i = 0; i < Ni - 1; ++i)
for (int i = 0; i < Ni; ++i)
{
Bx(i, j) -= FDTD_Const::C * dt / 2.0 * (Ez(i, j + 1) - Ez(i, j)) / dy;
By(i, j) += FDTD_Const::C * dt / 2.0 * (Ez(i + 1, j) - Ez(i, j)) / dx;
Bz(i, j) -= FDTD_Const::C * dt / 2.0 * ((Ey(i + 1, j) - Ey(i, j)) / dx - (Ex(i, j + 1) - Ex(i, j)) / dy);
}
}
periodical_boundary_conditions_shiftedB();

periodical_boundary_conditions_shiftedE();
#pragma omp parallel for
for (int j = 1; j < Nj; ++j)
for (int j = 0; j < Nj; ++j)
{
#pragma ivdep
for (int i = 1; i < Ni; ++i)
for (int i = 0; i < Ni; ++i)
{
Ex(i, j) += FDTD_Const::C * dt * (Bz(i, j) - Bz(i, j - 1)) / dy;
Ey(i, j) -= FDTD_Const::C * dt * (Bz(i, j) - Bz(i - 1, j)) / dx;
Expand All @@ -118,152 +121,15 @@ void FDTD::shifted_update_field(const double& time)
}

#pragma omp parallel for
for (int j = 0; j < Nj - 1; ++j)
for (int j = 0; j < Nj; ++j)
{
#pragma ivdep
for (int i = 0; i < Ni - 1; ++i)
for (int i = 0; i < Ni; ++i)
{
Bx(i, j) -= FDTD_Const::C * dt / 2.0 * (Ez(i, j + 1) - Ez(i, j)) / dy;
By(i, j) += FDTD_Const::C * dt / 2.0 * (Ez(i + 1, j) - Ez(i, j)) / dx;
Bz(i, j) -= FDTD_Const::C * dt / 2.0 * ((Ey(i + 1, j) - Ey(i, j)) / dx - (Ex(i, j + 1) - Ex(i, j)) / dy);
}
}
periodical_boundary_conditions_shiftedB();
}
}

void FDTD::periodical_boundary_conditions_shiftedE()
{
// Left boundary conditions
Ex(0, 0) += FDTD_Const::C * dt * (Bz(0, 0) - Bz(0, Nj - 1)) / dy;
Ey(0, 0) -= FDTD_Const::C * dt * (Bz(0, 0) - Bz(Ni - 1, 0)) / dx;
Ez(0, 0) += FDTD_Const::C * dt * ((By(0, 0) - By(Ni - 1, 0)) / dx - (Bx(0, 0) - Bx(0, Nj - 1)) / dy);
#pragma omp parallel for
for (int i = 1; i < Ni; ++i)
{
Ex(i, 0) += FDTD_Const::C * dt * (Bz(i, 0) - Bz(i, Nj - 1)) / dy;
Ey(i, 0) -= FDTD_Const::C * dt * (Bz(i, 0) - Bz(i - 1, 0)) / dx;
Ez(i, 0) += FDTD_Const::C * dt * ((By(i, 0) - By(i - 1, 0)) / dx - (Bx(i, 0) - Bx(i, Nj - 1)) / dy);
}

#pragma omp parallel for
for (int j = 1; j < Nj; ++j)
{
Ex(0, j) += FDTD_Const::C * dt * (Bz(0, j) - Bz(0, j - 1)) / dy;
Ey(0, j) -= FDTD_Const::C * dt * (Bz(0, j) - Bz(Ni - 1, j)) / dx;
Ez(0, j) += FDTD_Const::C * dt * ((By(0, j) - By(Ni - 1, j)) / dx - (Bx(0, j) - Bx(0, j - 1)) / dy);
}
}

void FDTD::periodical_boundary_conditions_shiftedB()
{
#pragma omp parallel for
for (int j = 0; j < Nj - 1; ++j)
{
Bx(Ni - 1, j) -= FDTD_Const::C * dt / 2.0 * (Ez(Ni - 1, j + 1) - Ez(Ni - 1, j)) / dy;
By(Ni - 1, j) += FDTD_Const::C * dt / 2.0 * (Ez(0, j) - Ez(Ni - 1, j)) / dx;
Bz(Ni - 1, j) -= FDTD_Const::C * dt / 2.0 * ((Ey(0, j) - Ey(Ni - 1, j)) / dx - (Ex(Ni - 1, j + 1) - Ex(Ni - 1, j)) / dy);
}

// Right boundary conditions
#pragma omp parallel for
for (int i = 0; i < Ni - 1; ++i)
{
Bx(i, Nj - 1) -= FDTD_Const::C * dt / 2.0 * (Ez(i, 0) - Ez(i, Nj - 1)) / dy;
By(i, Nj - 1) += FDTD_Const::C * dt / 2.0 * (Ez(i + 1, Nj - 1) - Ez(i, Nj - 1)) / dx;
Bz(i, Nj - 1) -= FDTD_Const::C * dt / 2.0 * ((Ey(i + 1, Nj - 1) - Ey(i, Nj - 1)) / dx - (Ex(i, 0) - Ex(i, Nj - 1)) / dy);
}
Bx(Ni - 1, Nj - 1) -= FDTD_Const::C * dt / 2.0 * (Ez(Ni - 1, 0) - Ez(Ni - 1, Nj - 1)) / dy;
By(Ni - 1, Nj - 1) += FDTD_Const::C * dt / 2.0 * (Ez(0, Nj - 1) - Ez(Ni - 1, Nj - 1)) / dx;
Bz(Ni - 1, Nj - 1) -= FDTD_Const::C * dt / 2.0 * ((Ey(0, Nj - 1) - Ey(Ni - 1, Nj - 1)) / dx - (Ex(Ni - 1, 0) - Ex(Ni - 1, Nj - 1)) / dy);
}

void FDTD::periodical_boundary_conditions_E()
{
// Left boundary conditions
Ex(0, 0) += FDTD_Const::C * dt * (Bz(0, 1) - Bz(0, Nj - 1)) / (2.0 * dy);
Ey(0, 0) -= FDTD_Const::C * dt * (Bz(1, 0) - Bz(Ni - 1, 0)) / (2.0 * dx);
Ez(0, 0) += FDTD_Const::C * dt * ((By(1, 0) - By(Ni - 1, 0)) / (2.0 * dx) - (Bx(0, 1) - Bx(0, Nj - 1)) / (2.0 * dy));
#pragma omp parallel for
for (int i = 1; i < Ni - 1; ++i)
{
Ex(i, 0) += FDTD_Const::C * dt * (Bz(i, 1) - Bz(i, Nj - 1)) / (2.0 * dy);
Ey(i, 0) -= FDTD_Const::C * dt * (Bz(i + 1, 0) - Bz(i - 1, 0)) / (2.0 * dx);
Ez(i, 0) += FDTD_Const::C * dt * ((By(i + 1, 0) - By(i - 1, 0)) / (2.0 * dx) - (Bx(i, 1) - Bx(i, Nj - 1)) / (2.0 * dy));
}
Ex(Ni - 1, 0) += FDTD_Const::C * dt * (Bz(Ni - 1, 1) - Bz(Ni - 1, Nj - 1)) / (2.0 * dy);
Ey(Ni - 1, 0) -= FDTD_Const::C * dt * (Bz(0, 0) - Bz(Ni - 2, 0)) / (2.0 * dx);
Ez(Ni - 1, 0) += FDTD_Const::C * dt * ((By(0, 0) - By(Ni - 2, 0)) / (2.0 * dx) - (Bx(Ni - 1, 1) - Bx(Ni - 1, Nj - 1)) / (2.0 * dy));

#pragma omp parallel for
for (int j = 1; j < Nj - 1; ++j)
{
Ex(0, j) += FDTD_Const::C * dt * (Bz(0, j + 1) - Bz(0, j - 1)) / (2.0 * dy);
Ey(0, j) -= FDTD_Const::C * dt * (Bz(1, j) - Bz(Ni - 1, j)) / (2.0 * dx);
Ez(0, j) += FDTD_Const::C * dt * ((By(1, j) - By(Ni - 1, j)) / (2.0 * dx) - (Bx(0, j + 1) - Bx(0, j - 1)) / (2.0 * dy));

Ex(Ni - 1, j) += FDTD_Const::C * dt * (Bz(Ni - 1, j + 1) - Bz(Ni - 1, j - 1)) / (2.0 * dy);
Ey(Ni - 1, j) -= FDTD_Const::C * dt * (Bz(0, j) - Bz(Ni - 2, j)) / (2.0 * dx);
Ez(Ni - 1, j) += FDTD_Const::C * dt * ((By(0, j) - By(Ni - 2, j)) / (2.0 * dx) - (Bx(Ni - 1, j + 1) - Bx(Ni - 1, j - 1)) / (2.0 * dy));
}

// Right boundary conditions
Ex(0, Nj - 1) += FDTD_Const::C * dt * (Bz(0, 0) - Bz(0, Nj - 2)) / (2.0 * dy);
Ey(0, Nj - 1) -= FDTD_Const::C * dt * (Bz(1, Nj - 1) - Bz(Ni - 1, Nj - 1)) / (2.0 * dx);
Ez(0, Nj - 1) += FDTD_Const::C * dt * ((By(1, Nj - 1) - By(Ni - 1, Nj - 1)) / (2.0 * dx) - (Bx(0, 0) - Bx(0, Nj - 2)) / (2.0 * dy));
#pragma omp parallel for
for (int i = 1; i < Ni - 1; ++i)
{
Ex(i, Nj - 1) += FDTD_Const::C * dt * (Bz(i, 0) - Bz(i, Nj - 2)) / (2.0 * dy);
Ey(i, Nj - 1) -= FDTD_Const::C * dt * (Bz(i + 1, Nj - 1) - Bz(i - 1, Nj - 1)) / (2.0 * dx);
Ez(i, Nj - 1) += FDTD_Const::C * dt * ((By(i + 1, Nj - 1) - By(i - 1, Nj - 1)) / (2.0 * dx) - (Bx(i, 0) - Bx(i, Nj - 2)) / (2.0 * dy));
}
Ex(Ni - 1, Nj - 1) += FDTD_Const::C * dt * (Bz(Ni - 1, 0) - Bz(Ni - 1, Nj - 2)) / (2.0 * dy);
Ey(Ni - 1, Nj - 1) -= FDTD_Const::C * dt * (Bz(0, Nj - 1) - Bz(Ni - 2, Nj - 1)) / (2.0 * dx);
Ez(Ni - 1, Nj - 1) += FDTD_Const::C * dt * ((By(0, Nj - 1) - By(Ni - 2, Nj - 1)) / (2.0 * dx) - (Bx(Ni - 1, 0) - Bx(Ni - 1, Nj - 2)) / (2.0 * dy));
}

void FDTD::periodical_boundary_conditions_B()
{
// Left boundary conditions
Bx(0, 0) -= FDTD_Const::C * dt * (Ez(0, 1) - Ez(0, Nj - 1)) / (2.0 * dy);
By(0, 0) += FDTD_Const::C * dt * (Ez(1, 0) - Ez(Ni - 1, Nj - 1)) / (2.0 * dx);
Bz(0, 0) -= FDTD_Const::C * dt * ((Ey(1, 0) - Ey(Ni - 1, 0)) / (2.0 * dx) - (Ex(0, 1) - Ex(0, Nj - 1)) / (2.0 * dy));
#pragma omp parallel for
for (int i = 1; i < Ni - 1; ++i)
{
Bx(i, 0) -= FDTD_Const::C * dt * (Ez(i, 1) - Ez(i, Nj - 1)) / (2.0 * dy);
By(i, 0) += FDTD_Const::C * dt * (Ez(i + 1, 0) - Ez(i - 1, 0)) / (2.0 * dx);
Bz(i, 0) -= FDTD_Const::C * dt * ((Ey(i + 1, 0) - Ey(i - 1, 0)) / (2.0 * dx) - (Ex(i, 1) - Ex(i, Nj - 1)) / (2.0 * dy));
}
Bx(Ni - 1, 0) -= FDTD_Const::C * dt * (Ez(Ni - 1, 1) - Ez(Ni - 1, Nj - 1)) / (2.0 * dy);
By(Ni - 1, 0) += FDTD_Const::C * dt * (Ez(0, 0) - Ez(Ni - 2, 0)) / (2.0 * dx);
Bz(Ni - 1, 0) -= FDTD_Const::C * dt * ((Ey(0, 0) - Ey(Ni - 2, 0)) / (2.0 * dx) - (Ex(Ni - 1, 1) - Ex(Ni - 1, Nj - 1)) / (2.0 * dy));

#pragma omp parallel for
for (int j = 1; j < Nj - 1; ++j)
{
Bx(0, j) -= FDTD_Const::C * dt * (Ez(0, j + 1) - Ez(0, j - 1)) / (2.0 * dy);
By(0, j) += FDTD_Const::C * dt * (Ez(1, j) - Ez(Ni - 1, j)) / (2.0 * dx);
Bz(0, j) -= FDTD_Const::C * dt * ((Ey(1, j) - Ey(Ni - 1, j)) / (2.0 * dx) - (Ex(0, j + 1) - Ex(0, j - 1)) / (2.0 * dy));

Bx(Ni - 1, j) -= FDTD_Const::C * dt * (Ez(Ni - 1, j + 1) - Ez(Ni - 1, j - 1)) / (2.0 * dy);
By(Ni - 1, j) += FDTD_Const::C * dt * (Ez(0, j) - Ez(Ni - 2, j)) / (2.0 * dx);
Bz(Ni - 1, j) -= FDTD_Const::C * dt * ((Ey(0, j) - Ey(Ni - 2, j)) / (2.0 * dx) - (Ex(Ni - 1, j + 1) - Ex(Ni - 1, j - 1)) / (2.0 * dy));
}

// Right boundary conditions
Bx(0, Nj - 1) -= FDTD_Const::C * dt * (Ez(0, 0) - Ez(0, Nj - 2)) / (2.0 * dy);
By(0, Nj - 1) += FDTD_Const::C * dt * (Ez(1, Nj - 1) - Ez(Ni - 1, Nj - 1)) / (2.0 * dx);
Bz(0, Nj - 1) -= FDTD_Const::C * dt * ((Ey(1, Nj - 1) - Ey(Ni - 1, Nj - 1)) / (2.0 * dx) - (Ex(0, 0) - Ex(0, Nj - 2)) / (2.0 * dy));
#pragma omp parallel for
for (int i = 1; i < Ni - 1; ++i)
{
Bx(i, Nj - 1) -= FDTD_Const::C * dt * (Ez(i, 0) - Ez(i, Nj - 2)) / (2.0 * dy);
By(i, Nj - 1) += FDTD_Const::C * dt * (Ez(i + 1, Nj - 1) - Ez(i - 1, Nj - 1)) / (2.0 * dx);
Bz(i, Nj - 1) -= FDTD_Const::C * dt * ((Ey(i + 1, Nj - 1) - Ey(i - 1, Nj - 1)) / (2.0 * dx) - (Ex(i, 0) - Ex(i, Nj - 2)) / (2.0 * dy));
}
Bx(Ni - 1, Nj - 1) -= FDTD_Const::C * dt * (Ez(Ni - 1, 0) - Ez(Ni - 1, Nj - 2)) / (2.0 * dy);
By(Ni - 1, Nj - 1) += FDTD_Const::C * dt * (Ez(0, Nj - 1) - Ez(Ni - 2, Nj - 1)) / (2.0 * dx);
Bz(Ni - 1, Nj - 1) -= FDTD_Const::C * dt * ((Ey(0, Nj - 1) - Ey(Ni - 2, Nj - 1)) / (2.0 * dx) - (Ex(Ni - 1, 0) - Ex(Ni - 1, Nj - 2)) / (2.0 * dy));
}
Loading

0 comments on commit 09ecacf

Please sign in to comment.