aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/BoundaryConditions/PML.cpp21
-rw-r--r--Source/BoundaryConditions/WarpXEvolvePML.cpp11
-rw-r--r--Source/BoundaryConditions/WarpX_PML_kernels.H76
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp60
-rw-r--r--Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H6
-rw-r--r--Source/FieldSolver/WarpXPushFieldsEM.cpp4
6 files changed, 121 insertions, 57 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp
index 38ea94364..8bc858615 100644
--- a/Source/BoundaryConditions/PML.cpp
+++ b/Source/BoundaryConditions/PML.cpp
@@ -495,12 +495,16 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
ngf = ngFFT;
}
+ // Allocate diagonal components (xx,yy,zz) only with divergence cleaning
+ const int ncomp = (do_dive_cleaning) ? 3 : 2;
+
pml_E_fp[0] = std::make_unique<MultiFab>(amrex::convert( ba,
- WarpX::GetInstance().getEfield_fp(0,0).ixType().toIntVect() ), dm, 3, nge );
+ WarpX::GetInstance().getEfield_fp(0,0).ixType().toIntVect() ), dm, ncomp, nge );
pml_E_fp[1] = std::make_unique<MultiFab>(amrex::convert( ba,
- WarpX::GetInstance().getEfield_fp(0,1).ixType().toIntVect() ), dm, 3, nge );
+ WarpX::GetInstance().getEfield_fp(0,1).ixType().toIntVect() ), dm, ncomp, nge );
pml_E_fp[2] = std::make_unique<MultiFab>(amrex::convert( ba,
- WarpX::GetInstance().getEfield_fp(0,2).ixType().toIntVect() ), dm, 3, nge );
+ WarpX::GetInstance().getEfield_fp(0,2).ixType().toIntVect() ), dm, ncomp, nge );
+
pml_B_fp[0] = std::make_unique<MultiFab>(amrex::convert( ba,
WarpX::GetInstance().getBfield_fp(0,0).ixType().toIntVect() ), dm, 2, ngb );
pml_B_fp[1] = std::make_unique<MultiFab>(amrex::convert( ba,
@@ -508,7 +512,6 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
pml_B_fp[2] = std::make_unique<MultiFab>(amrex::convert( ba,
WarpX::GetInstance().getBfield_fp(0,2).ixType().toIntVect() ), dm, 2, ngb );
-
pml_E_fp[0]->setVal(0.0);
pml_E_fp[1]->setVal(0.0);
pml_E_fp[2]->setVal(0.0);
@@ -522,6 +525,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
WarpX::GetInstance().getcurrent_fp(0,1).ixType().toIntVect() ), dm, 1, ngb );
pml_j_fp[2] = std::make_unique<MultiFab>(amrex::convert( ba,
WarpX::GetInstance().getcurrent_fp(0,2).ixType().toIntVect() ), dm, 1, ngb );
+
pml_j_fp[0]->setVal(0.0);
pml_j_fp[1]->setVal(0.0);
pml_j_fp[2]->setVal(0.0);
@@ -591,11 +595,12 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
DistributionMapping cdm{cba};
pml_E_cp[0] = std::make_unique<MultiFab>(amrex::convert( cba,
- WarpX::GetInstance().getEfield_cp(1,0).ixType().toIntVect() ), cdm, 3, nge );
+ WarpX::GetInstance().getEfield_cp(1,0).ixType().toIntVect() ), cdm, ncomp, nge );
pml_E_cp[1] = std::make_unique<MultiFab>(amrex::convert( cba,
- WarpX::GetInstance().getEfield_cp(1,1).ixType().toIntVect() ), cdm, 3, nge );
+ WarpX::GetInstance().getEfield_cp(1,1).ixType().toIntVect() ), cdm, ncomp, nge );
pml_E_cp[2] = std::make_unique<MultiFab>(amrex::convert( cba,
- WarpX::GetInstance().getEfield_cp(1,2).ixType().toIntVect() ), cdm, 3, nge );
+ WarpX::GetInstance().getEfield_cp(1,2).ixType().toIntVect() ), cdm, ncomp, nge );
+
pml_B_cp[0] = std::make_unique<MultiFab>(amrex::convert( cba,
WarpX::GetInstance().getBfield_cp(1,0).ixType().toIntVect() ), cdm, 2, ngb );
pml_B_cp[1] = std::make_unique<MultiFab>(amrex::convert( cba,
@@ -616,12 +621,14 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& /*grid_dm*/,
pml_F_cp->setVal(0.0);
}
+
pml_j_cp[0] = std::make_unique<MultiFab>(amrex::convert( cba,
WarpX::GetInstance().getcurrent_cp(1,0).ixType().toIntVect() ), cdm, 1, ngb );
pml_j_cp[1] = std::make_unique<MultiFab>(amrex::convert( cba,
WarpX::GetInstance().getcurrent_cp(1,1).ixType().toIntVect() ), cdm, 1, ngb );
pml_j_cp[2] = std::make_unique<MultiFab>(amrex::convert( cba,
WarpX::GetInstance().getcurrent_cp(1,2).ixType().toIntVect() ), cdm, 1, ngb );
+
pml_j_cp[0]->setVal(0.0);
pml_j_cp[1]->setVal(0.0);
pml_j_cp[2]->setVal(0.0);
diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp
index 31c8cbd1e..a72ab9afb 100644
--- a/Source/BoundaryConditions/WarpXEvolvePML.cpp
+++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp
@@ -43,6 +43,8 @@ WarpX::DampPML (int lev, PatchType patch_type)
WARPX_PROFILE("WarpX::DampPML()");
+ const bool dive_cleaning = WarpX::do_dive_cleaning;
+
if (pml[lev]->ok())
{
const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
@@ -110,15 +112,18 @@ WarpX::DampPML (int lev, PatchType patch_type)
amrex::ParallelFor(tex, tey, tez,
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
warpx_damp_pml_ex(i, j, k, pml_Exfab, Ex_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
- sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo);
+ sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
+ dive_cleaning);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
warpx_damp_pml_ey(i, j, k, pml_Eyfab, Ey_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
- sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo);
+ sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
+ dive_cleaning);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k) {
warpx_damp_pml_ez(i, j, k, pml_Ezfab, Ez_stag, sigma_fac_x, sigma_fac_y, sigma_fac_z,
- sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo);
+ sigma_star_fac_x, sigma_star_fac_y, sigma_star_fac_z, x_lo, y_lo, z_lo,
+ dive_cleaning);
});
amrex::ParallelFor(tbx, tby, tbz,
diff --git a/Source/BoundaryConditions/WarpX_PML_kernels.H b/Source/BoundaryConditions/WarpX_PML_kernels.H
index 7c29a6541..9e1b543e6 100644
--- a/Source/BoundaryConditions/WarpX_PML_kernels.H
+++ b/Source/BoundaryConditions/WarpX_PML_kernels.H
@@ -23,7 +23,8 @@ void warpx_damp_pml_ex (int i, int j, int k, Array4<Real> const& Ex,
const Real* const sigma_star_fac_x,
const Real* const sigma_star_fac_y,
const Real* const sigma_star_fac_z,
- int xlo, int ylo, int zlo)
+ int xlo, int ylo, int zlo,
+ const bool dive_cleaning)
{
#if (AMREX_SPACEDIM == 2)
@@ -33,11 +34,14 @@ void warpx_damp_pml_ex (int i, int j, int k, Array4<Real> const& Ex,
const int sx = Ex_stag[0];
const int sz = Ex_stag[1];
- // Exx
- if (sx == 0) {
- Ex(i,j,k,PMLComp::xx) *= sigma_star_fac_x[i-xlo];
- } else {
- Ex(i,j,k,PMLComp::xx) *= sigma_fac_x[i-xlo];
+ if (dive_cleaning)
+ {
+ // Exx
+ if (sx == 0) {
+ Ex(i,j,k,PMLComp::xx) *= sigma_star_fac_x[i-xlo];
+ } else {
+ Ex(i,j,k,PMLComp::xx) *= sigma_fac_x[i-xlo];
+ }
}
// Exz
@@ -54,11 +58,14 @@ void warpx_damp_pml_ex (int i, int j, int k, Array4<Real> const& Ex,
const int sy = Ex_stag[1];
const int sz = Ex_stag[2];
- // Exx
- if (sx == 0) {
- Ex(i,j,k,PMLComp::xx) *= sigma_star_fac_x[i-xlo];
- } else {
- Ex(i,j,k,PMLComp::xx) *= sigma_fac_x[i-xlo];
+ if (dive_cleaning)
+ {
+ // Exx
+ if (sx == 0) {
+ Ex(i,j,k,PMLComp::xx) *= sigma_star_fac_x[i-xlo];
+ } else {
+ Ex(i,j,k,PMLComp::xx) *= sigma_fac_x[i-xlo];
+ }
}
// Exy
@@ -87,11 +94,12 @@ void warpx_damp_pml_ey (int i, int j, int k, Array4<Real> const& Ey,
const Real* const sigma_star_fac_x,
const Real* const sigma_star_fac_y,
const Real* const sigma_star_fac_z,
- int xlo, int ylo, int zlo)
+ int xlo, int ylo, int zlo,
+ const bool dive_cleaning)
{
#if (AMREX_SPACEDIM == 2)
- amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo);
+ amrex::ignore_unused(sigma_fac_y, sigma_star_fac_y, ylo, dive_cleaning);
// sx = 0 means that Ey is staggered in x, while sx = 1 means that Ey is nodal in x (same for z)
const int sx = Ey_stag[0];
@@ -125,11 +133,14 @@ void warpx_damp_pml_ey (int i, int j, int k, Array4<Real> const& Ey,
Ey(i,j,k,PMLComp::yx) *= sigma_fac_x[i-xlo];
}
- // Eyy
- if (sy == 0) {
- Ey(i,j,k,PMLComp::yy) *= sigma_star_fac_y[j-ylo];
- } else {
- Ey(i,j,k,PMLComp::yy) *= sigma_fac_y[j-ylo];
+ if (dive_cleaning)
+ {
+ // Eyy
+ if (sy == 0) {
+ Ey(i,j,k,PMLComp::yy) *= sigma_star_fac_y[j-ylo];
+ } else {
+ Ey(i,j,k,PMLComp::yy) *= sigma_fac_y[j-ylo];
+ }
}
// Eyz
@@ -151,7 +162,8 @@ void warpx_damp_pml_ez (int i, int j, int k, Array4<Real> const& Ez,
const Real* const sigma_star_fac_x,
const Real* const sigma_star_fac_y,
const Real* const sigma_star_fac_z,
- int xlo, int ylo, int zlo)
+ int xlo, int ylo, int zlo,
+ const bool dive_cleaning)
{
#if (AMREX_SPACEDIM == 2)
@@ -168,11 +180,14 @@ void warpx_damp_pml_ez (int i, int j, int k, Array4<Real> const& Ez,
Ez(i,j,k,PMLComp::zx) *= sigma_fac_x[i-xlo];
}
- // Ezz
- if (sz == 0) {
- Ez(i,j,k,PMLComp::zz) *= sigma_star_fac_z[j-zlo];
- } else {
- Ez(i,j,k,PMLComp::zz) *= sigma_fac_z[j-zlo];
+ if (dive_cleaning)
+ {
+ // Ezz
+ if (sz == 0) {
+ Ez(i,j,k,PMLComp::zz) *= sigma_star_fac_z[j-zlo];
+ } else {
+ Ez(i,j,k,PMLComp::zz) *= sigma_fac_z[j-zlo];
+ }
}
#elif (AMREX_SPACEDIM == 3)
@@ -196,11 +211,14 @@ void warpx_damp_pml_ez (int i, int j, int k, Array4<Real> const& Ez,
Ez(i,j,k,PMLComp::zy) *= sigma_fac_y[j-ylo];
}
- // Ezz
- if (sz == 0) {
- Ez(i,j,k,PMLComp::zz) *= sigma_star_fac_z[k-zlo];
- } else {
- Ez(i,j,k,PMLComp::zz) *= sigma_fac_z[k-zlo];
+ if (dive_cleaning)
+ {
+ // Ezz
+ if (sz == 0) {
+ Ez(i,j,k,PMLComp::zz) *= sigma_star_fac_z[k-zlo];
+ } else {
+ Ez(i,j,k,PMLComp::zz) *= sigma_fac_z[k-zlo];
+ }
}
#endif
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
index 8346eb067..6fc157653 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveBPML.cpp
@@ -26,7 +26,8 @@ using namespace amrex;
void FiniteDifferenceSolver::EvolveBPML (
std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > const Efield,
- amrex::Real const dt ) {
+ amrex::Real const dt,
+ const bool dive_cleaning) {
// Select algorithm (The choice of algorithm is a runtime option,
// but we compile code for each algorithm, using templates)
@@ -36,15 +37,15 @@ void FiniteDifferenceSolver::EvolveBPML (
#else
if (m_do_nodal) {
- EvolveBPMLCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, dt );
+ EvolveBPMLCartesian <CartesianNodalAlgorithm> (Bfield, Efield, dt, dive_cleaning);
} else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
- EvolveBPMLCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, dt );
+ EvolveBPMLCartesian <CartesianYeeAlgorithm> (Bfield, Efield, dt, dive_cleaning);
} else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
- EvolveBPMLCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, dt );
+ EvolveBPMLCartesian <CartesianCKCAlgorithm> (Bfield, Efield, dt, dive_cleaning);
} else {
amrex::Abort("EvolveBPML: Unknown algorithm");
@@ -59,7 +60,8 @@ template<typename T_Algo>
void FiniteDifferenceSolver::EvolveBPMLCartesian (
std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > const Efield,
- amrex::Real const dt ) {
+ amrex::Real const dt,
+ const bool dive_cleaning) {
// Loop through the grids, and over the tiles within each grid
#ifdef _OPENMP
@@ -92,36 +94,66 @@ void FiniteDifferenceSolver::EvolveBPMLCartesian (
amrex::ParallelFor(tbx, tby, tbz,
[=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+ amrex::Real UpwardDz_Ey_yy = 0._rt;
+ amrex::Real UpwardDy_Ez_zz = 0._rt;
+ if (dive_cleaning)
+ {
+ UpwardDz_Ey_yy = T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yy);
+ UpwardDy_Ez_zz = T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz);
+ }
+
Bx(i, j, k, PMLComp::xz) += dt * (
T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yx)
- + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yy)
- + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz) );
+ + T_Algo::UpwardDz(Ey, coefs_z, n_coefs_z, i, j, k, PMLComp::yz)
+ + UpwardDz_Ey_yy);
+
Bx(i, j, k, PMLComp::xy) -= dt * (
T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zx)
+ T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zy)
- + T_Algo::UpwardDy(Ez, coefs_y, n_coefs_y, i, j, k, PMLComp::zz) );
+ + UpwardDy_Ez_zz);
},
[=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+ amrex::Real UpwardDx_Ez_zz = 0._rt;
+ amrex::Real UpwardDz_Ex_xx = 0._rt;
+ if (dive_cleaning)
+ {
+ UpwardDx_Ez_zz = T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz);
+ UpwardDz_Ex_xx = T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xx);
+ }
+
By(i, j, k, PMLComp::yx) += dt * (
T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zx)
+ T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zy)
- + T_Algo::UpwardDx(Ez, coefs_x, n_coefs_x, i, j, k, PMLComp::zz) );
+ + UpwardDx_Ez_zz);
+
By(i, j, k, PMLComp::yz) -= dt * (
- T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xx)
+ UpwardDz_Ex_xx
+ T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xy)
- + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xz) );
+ + T_Algo::UpwardDz(Ex, coefs_z, n_coefs_z, i, j, k, PMLComp::xz));
},
[=] AMREX_GPU_DEVICE (int i, int j, int k){
+
+ amrex::Real UpwardDy_Ex_xx = 0._rt;
+ amrex::Real UpwardDx_Ey_yy = 0._rt;
+ if (dive_cleaning)
+ {
+ UpwardDy_Ex_xx = T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xx);
+ UpwardDx_Ey_yy = T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yy);
+ }
+
Bz(i, j, k, PMLComp::zy) += dt * (
- T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xx)
+ UpwardDy_Ex_xx
+ T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xy)
+ T_Algo::UpwardDy(Ex, coefs_y, n_coefs_y, i, j, k, PMLComp::xz) );
+
Bz(i, j, k, PMLComp::zx) -= dt * (
T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yx)
- + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yy)
- + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz) );
+ + T_Algo::UpwardDx(Ey, coefs_x, n_coefs_x, i, j, k, PMLComp::yz)
+ + UpwardDx_Ey_yy);
}
);
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index a311bb2bc..bb96dee45 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -75,7 +75,8 @@ class FiniteDifferenceSolver
void EvolveBPML ( std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > const Efield,
- amrex::Real const dt );
+ amrex::Real const dt,
+ const bool dive_cleaning);
void EvolveEPML ( std::array< amrex::MultiFab*, 3 > Efield,
std::array< amrex::MultiFab*, 3 > const Bfield,
@@ -176,7 +177,8 @@ class FiniteDifferenceSolver
void EvolveBPMLCartesian (
std::array< amrex::MultiFab*, 3 > Bfield,
std::array< amrex::MultiFab*, 3 > const Efield,
- amrex::Real const dt );
+ amrex::Real const dt,
+ const bool dive_cleaning);
template< typename T_Algo >
void EvolveEPMLCartesian (
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 45b7b3b99..63d519eee 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -192,10 +192,10 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt)
if (do_pml && pml[lev]->ok()) {
if (patch_type == PatchType::fine) {
m_fdtd_solver_fp[lev]->EvolveBPML(
- pml[lev]->GetB_fp(), pml[lev]->GetE_fp(), a_dt );
+ pml[lev]->GetB_fp(), pml[lev]->GetE_fp(), a_dt, WarpX::do_dive_cleaning);
} else {
m_fdtd_solver_cp[lev]->EvolveBPML(
- pml[lev]->GetB_cp(), pml[lev]->GetE_cp(), a_dt );
+ pml[lev]->GetB_cp(), pml[lev]->GetE_cp(), a_dt, WarpX::do_dive_cleaning);
}
}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-12 10:57:22 +0000