Skip to content

Commit

Permalink
update pml
Browse files Browse the repository at this point in the history
  • Loading branch information
@Amazingkivas
Amazingkivas committed May 4, 2024
1 parent a2da572 commit 34db933
Show file tree
Hide file tree
Showing 6 changed files with 83 additions and 73 deletions.
Binary file added PlotScript/Animations/animation_Bx.gif
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified PlotScript/Animations/animation_Ex.gif
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified PlotScript/src/Release/sample.exe
View file
Binary file not shown.
6 changes: 2 additions & 4 deletions PlotScript/visualization.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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'}

Expand All @@ -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)
Expand Down Expand Up @@ -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)
4 changes: 3 additions & 1 deletion samples/sample.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>(n) / 2.0 * d;

Parameters params
{
n,
Expand Down Expand Up @@ -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

Expand Down
146 changes: 78 additions & 68 deletions src/FDTD.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ void FDTD::set_sigma_x(int bounds_i[2], int bounds_j[2], int bounds_k[2],
pow((static_cast<double>(dist(i, j, k))) /
static_cast<double>(pml_size_i), FDTDconst::N);
BsigmaX(i, j, k) = SGm *
pow((static_cast<double>(dist(i, j, k)) + 0.5) /
pow((static_cast<double>(dist(i, j, k))) /
static_cast<double>(pml_size_i), FDTDconst::N);
}
}
Expand All @@ -78,7 +78,7 @@ void FDTD::set_sigma_y(int bounds_i[2], int bounds_j[2], int bounds_k[2],
pow((static_cast<double>(dist(i, j, k))) /
static_cast<double>(pml_size_j), FDTDconst::N);
BsigmaY(i, j, k) = SGm *
pow((static_cast<double>(dist(i, j, k)) + 0.5) /
pow((static_cast<double>(dist(i, j, k))) /
static_cast<double>(pml_size_j), FDTDconst::N);
}
}
Expand All @@ -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<double>(dist(i, j, k)))
/ static_cast<double>(pml_size_k), FDTDconst::N);
pow((static_cast<double>(dist(i, j, k))) /
static_cast<double>(pml_size_k), FDTDconst::N);
BsigmaZ(i, j, k) = SGm *
pow((static_cast<double>(dist(i, j, k)) + 0.5)
/ static_cast<double>(pml_size_k), FDTDconst::N);
pow((static_cast<double>(dist(i, j, k))) /
static_cast<double>(pml_size_k), FDTDconst::N);
}
}
}
Expand All @@ -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 ||
Expand Down Expand Up @@ -234,40 +228,45 @@ 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++)
{
for (int j = bounds_j[0]; j < bounds_j[1]; j++)
{
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);
Expand All @@ -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++)
Expand All @@ -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);
Expand Down Expand Up @@ -519,17 +523,23 @@ std::vector<std::vector<Field>> FDTD::update_fields(const int time)
{
max_val_2 = calc_max_value(Ex);
}
#ifdef __OUTPUT_SIGMA__
if (t == 0) break;
std::vector<Field> new_iteration{ EsigmaX, EsigmaY, EsigmaZ,
BsigmaX, BsigmaY, BsigmaZ };
#endif // __OUTPUT_SIGMA__

#endif //__DEBUG_PML__

#ifndef __OUTPUT_SIGMA__
std::vector<Field> new_iteration{ Ex, Ey, Ez, Bx, By, Bz };
#endif // __OUTPUT_SIGMA__

#ifdef __OUTPUT_SIGMA__

EsigmaX(parameters.Ni / 2, 1, 1) = 100;

std::vector<Field> 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);
}

Expand Down

0 comments on commit 34db933

Please sign in to comment.