Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updating PML #9

Merged
merged 2 commits into from
May 4, 2024
Merged

Updating PML #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
34db933
update pml
Amazingkivas May 4, 2024
284ff11
minor fix
Amazingkivas May 4, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file added PlotScript/Animations/animation_Bx.gif
View file
Open in desktop
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
Open in desktop
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
Open in desktop
Binary file not shown.
4 changes: 1 addition & 3 deletions PlotScript/visualization.py
View file
Open in desktop
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 @@ -52,4 +51,3 @@ def execute_cpp(grid_size, iters_num, single_iteration_flag=True):
if __name__ == '__main__':
execute_cpp(75, 410, False)
get_animation("Ex")
# get_heatmap("Ex", 86)
4 changes: 3 additions & 1 deletion samples/sample.cpp
View file
Open in desktop
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
Open in desktop
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