aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Source/WarpX.H15
-rw-r--r--Source/WarpXComm.cpp137
-rw-r--r--Source/WarpXEvolve.cpp630
-rw-r--r--Source/WarpXPML.H10
-rw-r--r--Source/WarpXPML.cpp138
5 files changed, 524 insertions, 406 deletions
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 19b3d74ad..c0ab7495f 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -34,6 +34,12 @@ enum struct DtType : int
SecondHalf
};
+enum struct PatchType : int
+{
+ fine,
+ coarse
+};
+
class WarpX
: public amrex::AmrCore
{
@@ -230,6 +236,15 @@ private:
/// Advance the simulation by numsteps steps, electromagnetic case.
///
void EvolveEM(int numsteps);
+
+ void FillBoundaryB (int lev, PatchType patch_type);
+ void FillBoundaryE (int lev, PatchType patch_type);
+ void FillBoundaryF (int lev, PatchType patch_type);
+
+ void EvolveB (int lev, PatchType patch_type, amrex::Real dt);
+ void EvolveE (int lev, PatchType patch_type, amrex::Real dt);
+ void EvolveF (int lev, PatchType patch_type, amrex::Real dt, DtType dt_type);
+ void DampPML (int lev, PatchType patch_type);
#ifdef WARPX_DO_ELECTROSTATIC
///
diff --git a/Source/WarpXComm.cpp b/Source/WarpXComm.cpp
index 23ce0ce0b..46d3a0af7 100644
--- a/Source/WarpXComm.cpp
+++ b/Source/WarpXComm.cpp
@@ -238,67 +238,122 @@ WarpX::FillBoundaryF ()
void
WarpX::FillBoundaryE(int lev)
{
- const auto& period = Geom(lev).periodicity();
-
- if (do_pml && pml[lev]->ok()) {
- ExchangeWithPmlE(lev);
- pml[lev]->FillBoundaryE();
- }
-
- (*Efield_fp[lev][0]).FillBoundary( period );
- (*Efield_fp[lev][1]).FillBoundary( period );
- (*Efield_fp[lev][2]).FillBoundary( period );
+ FillBoundaryE(lev, PatchType::fine);
+ if (lev > 0) FillBoundaryE(lev, PatchType::coarse);
+}
- if (lev > 0)
+void
+WarpX::FillBoundaryE (int lev, PatchType patch_type)
+{
+ if (patch_type == PatchType::fine)
{
- const auto& cperiod = Geom(lev-1).periodicity();
- (*Efield_cp[lev][0]).FillBoundary(cperiod);
- (*Efield_cp[lev][1]).FillBoundary(cperiod);
- (*Efield_cp[lev][2]).FillBoundary(cperiod);
+ if (do_pml && pml[lev]->ok())
+ {
+ pml[lev]->ExchangeE(patch_type,
+ { Efield_fp[lev][0].get(),
+ Efield_fp[lev][1].get(),
+ Efield_fp[lev][2].get() });
+ pml[lev]->FillBoundaryE(patch_type);
+ }
+
+ const auto& period = Geom(lev).periodicity();
+ (*Efield_fp[lev][0]).FillBoundary(period);
+ (*Efield_fp[lev][1]).FillBoundary(period);
+ (*Efield_fp[lev][2]).FillBoundary(period);
}
+ else if (patch_type == PatchType::coarse)
+ {
+ if (do_pml && pml[lev]->ok())
+ {
+ pml[lev]->ExchangeE(patch_type,
+ { Efield_cp[lev][0].get(),
+ Efield_cp[lev][1].get(),
+ Efield_cp[lev][2].get() });
+ pml[lev]->FillBoundaryE(patch_type);
+ }
+
+ const auto& cperiod = Geom(lev-1).periodicity();
+ (*Efield_cp[lev][0]).FillBoundary(cperiod);
+ (*Efield_cp[lev][1]).FillBoundary(cperiod);
+ (*Efield_cp[lev][2]).FillBoundary(cperiod);
+ }
}
void
-WarpX::FillBoundaryB(int lev)
+WarpX::FillBoundaryB (int lev)
{
- const auto& period = Geom(lev).periodicity();
+ FillBoundaryB(lev, PatchType::fine);
+ if (lev > 0) FillBoundaryB(lev, PatchType::coarse);
+}
- if (do_pml && pml[lev]->ok())
+void
+WarpX::FillBoundaryB (int lev, PatchType patch_type)
+{
+ if (patch_type == PatchType::fine)
{
- ExchangeWithPmlB(lev);
- pml[lev]->FillBoundaryB();
- }
-
- (*Bfield_fp[lev][0]).FillBoundary(period);
- (*Bfield_fp[lev][1]).FillBoundary(period);
- (*Bfield_fp[lev][2]).FillBoundary(period);
+ if (do_pml && pml[lev]->ok())
+ {
+ pml[lev]->ExchangeB(patch_type,
+ { Bfield_fp[lev][0].get(),
+ Bfield_fp[lev][1].get(),
+ Bfield_fp[lev][2].get() });
+ pml[lev]->FillBoundaryB(patch_type);
+ }
- if (lev > 0)
+ const auto& period = Geom(lev).periodicity();
+ (*Bfield_fp[lev][0]).FillBoundary(period);
+ (*Bfield_fp[lev][1]).FillBoundary(period);
+ (*Bfield_fp[lev][2]).FillBoundary(period);
+ }
+ else if (patch_type == PatchType::coarse)
{
- const auto& cperiod = Geom(lev-1).periodicity();
- (*Bfield_cp[lev][0]).FillBoundary(cperiod);
- (*Bfield_cp[lev][1]).FillBoundary(cperiod);
- (*Bfield_cp[lev][2]).FillBoundary(cperiod);
+ if (do_pml && pml[lev]->ok())
+ {
+ pml[lev]->ExchangeB(patch_type,
+ { Bfield_cp[lev][0].get(),
+ Bfield_cp[lev][1].get(),
+ Bfield_cp[lev][2].get() });
+ pml[lev]->FillBoundaryB(patch_type);
+ }
+
+ const auto& cperiod = Geom(lev-1).periodicity();
+ (*Bfield_cp[lev][0]).FillBoundary(cperiod);
+ (*Bfield_cp[lev][1]).FillBoundary(cperiod);
+ (*Bfield_cp[lev][2]).FillBoundary(cperiod);
}
}
void
-WarpX::FillBoundaryF(int lev)
+WarpX::FillBoundaryF (int lev)
{
- const auto& period = Geom(lev).periodicity();
+ FillBoundaryF(lev, PatchType::fine);
+ if (lev > 0) FillBoundaryF(lev, PatchType::coarse);
+}
- if (do_pml && pml[lev]->ok())
+void
+WarpX::FillBoundaryF (int lev, PatchType patch_type)
+{
+ if (patch_type == PatchType::fine && F_fp[lev])
{
- ExchangeWithPmlF(lev);
- pml[lev]->FillBoundaryF();
- }
-
- if (F_fp[lev]) F_fp[lev]->FillBoundary(period);
+ if (do_pml && pml[lev]->ok())
+ {
+ pml[lev]->ExchangeF(patch_type, F_fp[lev].get());
+ pml[lev]->FillBoundaryF(patch_type);
+ }
- if (lev > 0)
+ const auto& period = Geom(lev).periodicity();
+ F_fp[lev]->FillBoundary(period);
+ }
+ else if (patch_type == PatchType::coarse && F_cp[lev])
{
- const auto& cperiod = Geom(lev-1).periodicity();
- if (F_cp[lev]) F_cp[lev]->FillBoundary(cperiod);
+ if (do_pml && pml[lev]->ok())
+ {
+ pml[lev]->ExchangeF(patch_type, F_cp[lev].get());
+ pml[lev]->FillBoundaryF(patch_type);
+ }
+
+ const auto& cperiod = Geom(lev-1).periodicity();
+ F_cp[lev]->FillBoundary(cperiod);
}
}
diff --git a/Source/WarpXEvolve.cpp b/Source/WarpXEvolve.cpp
index be5152fc2..f1a0faa15 100644
--- a/Source/WarpXEvolve.cpp
+++ b/Source/WarpXEvolve.cpp
@@ -262,107 +262,102 @@ void
WarpX::EvolveB (int lev, Real dt)
{
BL_PROFILE("WarpX::EvolveB()");
+ EvolveB(lev, PatchType::fine, dt);
+ if (lev > 0) {
+ EvolveB(lev, PatchType::coarse, dt);
+ }
+}
- // Parameters of the solver: order and mesh spacing
- const int norder = 2;
-
- int npatches = (lev == 0) ? 1 : 2;
-
- for (int ipatch = 0; ipatch < npatches; ++ipatch)
+void
+WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt)
+{
+ const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+ const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
+ const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
+
+ MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz;
+ if (patch_type == PatchType::fine)
{
- int patch_level = (ipatch == 0) ? lev : lev-1;
- const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
- const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
-
- MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz;
- if (ipatch == 0)
- {
- Ex = Efield_fp[lev][0].get();
- Ey = Efield_fp[lev][1].get();
- Ez = Efield_fp[lev][2].get();
- Bx = Bfield_fp[lev][0].get();
- By = Bfield_fp[lev][1].get();
- Bz = Bfield_fp[lev][2].get();
- }
- else
- {
- Ex = Efield_cp[lev][0].get();
- Ey = Efield_cp[lev][1].get();
- Ez = Efield_cp[lev][2].get();
- Bx = Bfield_cp[lev][0].get();
- By = Bfield_cp[lev][1].get();
- Bz = Bfield_cp[lev][2].get();
- }
+ Ex = Efield_fp[lev][0].get();
+ Ey = Efield_fp[lev][1].get();
+ Ez = Efield_fp[lev][2].get();
+ Bx = Bfield_fp[lev][0].get();
+ By = Bfield_fp[lev][1].get();
+ Bz = Bfield_fp[lev][2].get();
+ }
+ else
+ {
+ Ex = Efield_cp[lev][0].get();
+ Ey = Efield_cp[lev][1].get();
+ Ez = Efield_cp[lev][2].get();
+ Bx = Bfield_cp[lev][0].get();
+ By = Bfield_cp[lev][1].get();
+ Bz = Bfield_cp[lev][2].get();
+ }
- MultiFab* cost = costs[lev].get();
- const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector();
+ MultiFab* cost = costs[lev].get();
+ const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector();
- // Loop through the grids, and over the tiles within each grid
+ // Loop through the grids, and over the tiles within each grid
#ifdef _OPENMP
#pragma omp parallel
#endif
- for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi )
- {
- Real wt = amrex::second();
-
- const Box& tbx = mfi.tilebox(Bx_nodal_flag);
- const Box& tby = mfi.tilebox(By_nodal_flag);
- const Box& tbz = mfi.tilebox(Bz_nodal_flag);
-
- // Call picsar routine for each tile
- warpx_push_bvec(
- tbx.loVect(), tbx.hiVect(),
- tby.loVect(), tby.hiVect(),
- tbz.loVect(), tbz.hiVect(),
- BL_TO_FORTRAN_3D((*Ex)[mfi]),
- BL_TO_FORTRAN_3D((*Ey)[mfi]),
- BL_TO_FORTRAN_3D((*Ez)[mfi]),
- BL_TO_FORTRAN_3D((*Bx)[mfi]),
- BL_TO_FORTRAN_3D((*By)[mfi]),
- BL_TO_FORTRAN_3D((*Bz)[mfi]),
- &dtsdx[0], &dtsdx[1], &dtsdx[2],
- &WarpX::maxwell_fdtd_solver_id);
-
- if (cost) {
- Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
- if (ipatch == 1) cbx.refine(rr);
- wt = (amrex::second() - wt) / cbx.d_numPts();
- (*cost)[mfi].plus(wt, cbx);
- }
- }
+ for ( MFIter mfi(*Bx,true); mfi.isValid(); ++mfi )
+ {
+ Real wt = amrex::second();
+
+ const Box& tbx = mfi.tilebox(Bx_nodal_flag);
+ const Box& tby = mfi.tilebox(By_nodal_flag);
+ const Box& tbz = mfi.tilebox(Bz_nodal_flag);
+
+ // Call picsar routine for each tile
+ warpx_push_bvec(
+ tbx.loVect(), tbx.hiVect(),
+ tby.loVect(), tby.hiVect(),
+ tbz.loVect(), tbz.hiVect(),
+ BL_TO_FORTRAN_3D((*Ex)[mfi]),
+ BL_TO_FORTRAN_3D((*Ey)[mfi]),
+ BL_TO_FORTRAN_3D((*Ez)[mfi]),
+ BL_TO_FORTRAN_3D((*Bx)[mfi]),
+ BL_TO_FORTRAN_3D((*By)[mfi]),
+ BL_TO_FORTRAN_3D((*Bz)[mfi]),
+ &dtsdx[0], &dtsdx[1], &dtsdx[2],
+ &WarpX::maxwell_fdtd_solver_id);
+
+ if (cost) {
+ Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
+ if (patch_type == PatchType::coarse) cbx.refine(rr);
+ wt = (amrex::second() - wt) / cbx.d_numPts();
+ (*cost)[mfi].plus(wt, cbx);
+ }
}
- if (do_pml && pml[lev]->ok())
+ if (do_pml && pml[lev]->ok())
{
- for (int ipatch = 0; ipatch < npatches; ++ipatch)
- {
- const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
- const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
- int patch_level = (ipatch == 0) ? lev : lev-1;
- const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
- const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
+ const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
+ const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+
#ifdef _OPENMP
#pragma omp parallel
#endif
- for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi )
- {
- const Box& tbx = mfi.tilebox(Bx_nodal_flag);
- const Box& tby = mfi.tilebox(By_nodal_flag);
- const Box& tbz = mfi.tilebox(Bz_nodal_flag);
-
- WRPX_PUSH_PML_BVEC(
- tbx.loVect(), tbx.hiVect(),
- tby.loVect(), tby.hiVect(),
- tbz.loVect(), tbz.hiVect(),
- BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
- &dtsdx[0], &dtsdx[1], &dtsdx[2],
- &WarpX::maxwell_fdtd_solver_id);
- }
+ for ( MFIter mfi(*pml_B[0],true); mfi.isValid(); ++mfi )
+ {
+ const Box& tbx = mfi.tilebox(Bx_nodal_flag);
+ const Box& tby = mfi.tilebox(By_nodal_flag);
+ const Box& tbz = mfi.tilebox(Bz_nodal_flag);
+
+ WRPX_PUSH_PML_BVEC(
+ tbx.loVect(), tbx.hiVect(),
+ tby.loVect(), tby.hiVect(),
+ tbz.loVect(), tbz.hiVect(),
+ BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
+ &dtsdx[0], &dtsdx[1], &dtsdx[2],
+ &WarpX::maxwell_fdtd_solver_id);
}
}
}
@@ -370,236 +365,235 @@ WarpX::EvolveB (int lev, Real dt)
void
WarpX::EvolveE (Real dt)
{
- for (int lev = 0; lev <= finest_level; ++lev) {
- EvolveE(lev, dt);
- }
+ for (int lev = 0; lev <= finest_level; ++lev) {
+ EvolveE(lev, dt);
+ }
}
void
WarpX::EvolveE (int lev, Real dt)
{
- BL_PROFILE("WarpX::EvolveE()");
+ BL_PROFILE("WarpX::EvolveE()");
+ EvolveE(lev, PatchType::fine, dt);
+ if (lev > 0) {
+ EvolveE(lev, PatchType::coarse, dt);
+ }
+}
- // Parameters of the solver: order and mesh spacing
- const int norder = 2;
- static constexpr Real c2 = PhysConst::c*PhysConst::c;
- const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt;
- const Real c2dt = (PhysConst::c*PhysConst::c) * dt;
+void
+WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt)
+{
+ // Parameters of the solver: order and mesh spacing
+ static constexpr Real c2 = PhysConst::c*PhysConst::c;
+ const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt;
+ const Real c2dt = (PhysConst::c*PhysConst::c) * dt;
- int npatches = (lev == 0) ? 1 : 2;
+ int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+ const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
+ const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
- for (int ipatch = 0; ipatch < npatches; ++ipatch)
+ MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F;
+ if (patch_type == PatchType::fine)
{
- int patch_level = (ipatch == 0) ? lev : lev-1;
- const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
- const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
-
- MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz, *jx, *jy, *jz, *F;
- if (ipatch == 0)
- {
- Ex = Efield_fp[lev][0].get();
- Ey = Efield_fp[lev][1].get();
- Ez = Efield_fp[lev][2].get();
- Bx = Bfield_fp[lev][0].get();
- By = Bfield_fp[lev][1].get();
- Bz = Bfield_fp[lev][2].get();
- jx = current_fp[lev][0].get();
- jy = current_fp[lev][1].get();
- jz = current_fp[lev][2].get();
- F = F_fp[lev].get();
- }
- else
- {
- Ex = Efield_cp[lev][0].get();
- Ey = Efield_cp[lev][1].get();
- Ez = Efield_cp[lev][2].get();
- Bx = Bfield_cp[lev][0].get();
- By = Bfield_cp[lev][1].get();
- Bz = Bfield_cp[lev][2].get();
- jx = current_cp[lev][0].get();
- jy = current_cp[lev][1].get();
- jz = current_cp[lev][2].get();
- F = F_cp[lev].get();
- }
+ Ex = Efield_fp[lev][0].get();
+ Ey = Efield_fp[lev][1].get();
+ Ez = Efield_fp[lev][2].get();
+ Bx = Bfield_fp[lev][0].get();
+ By = Bfield_fp[lev][1].get();
+ Bz = Bfield_fp[lev][2].get();
+ jx = current_fp[lev][0].get();
+ jy = current_fp[lev][1].get();
+ jz = current_fp[lev][2].get();
+ F = F_fp[lev].get();
+ }
+ else if (patch_type == PatchType::coarse)
+ {
+ Ex = Efield_cp[lev][0].get();
+ Ey = Efield_cp[lev][1].get();
+ Ez = Efield_cp[lev][2].get();
+ Bx = Bfield_cp[lev][0].get();
+ By = Bfield_cp[lev][1].get();
+ Bz = Bfield_cp[lev][2].get();
+ jx = current_cp[lev][0].get();
+ jy = current_cp[lev][1].get();
+ jz = current_cp[lev][2].get();
+ F = F_cp[lev].get();
+ }
- MultiFab* cost = costs[lev].get();
- const IntVect& rr = (lev < finestLevel()) ? refRatio(lev) : IntVect::TheUnitVector();
+ MultiFab* cost = costs[lev].get();
+ const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector();
- // Loop through the grids, and over the tiles within each grid
+ // Loop through the grids, and over the tiles within each grid
#ifdef _OPENMP
#pragma omp parallel
#endif
- for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi )
- {
- Real wt = amrex::second();
-
- const Box& tex = mfi.tilebox(Ex_nodal_flag);
- const Box& tey = mfi.tilebox(Ey_nodal_flag);
- const Box& tez = mfi.tilebox(Ez_nodal_flag);
-
- // Call picsar routine for each tile
- warpx_push_evec(
- tex.loVect(), tex.hiVect(),
- tey.loVect(), tey.hiVect(),
- tez.loVect(), tez.hiVect(),
- BL_TO_FORTRAN_3D((*Ex)[mfi]),
- BL_TO_FORTRAN_3D((*Ey)[mfi]),
- BL_TO_FORTRAN_3D((*Ez)[mfi]),
- BL_TO_FORTRAN_3D((*Bx)[mfi]),
- BL_TO_FORTRAN_3D((*By)[mfi]),
- BL_TO_FORTRAN_3D((*Bz)[mfi]),
- BL_TO_FORTRAN_3D((*jx)[mfi]),
- BL_TO_FORTRAN_3D((*jy)[mfi]),
- BL_TO_FORTRAN_3D((*jz)[mfi]),
- &mu_c2_dt,
- &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
-
- if (F) {
-
- // Call picsar routine for each tile
- warpx_push_evec_f(
- tex.loVect(), tex.hiVect(),
- tey.loVect(), tey.hiVect(),
- tez.loVect(), tez.hiVect(),
- BL_TO_FORTRAN_3D((*Ex)[mfi]),
- BL_TO_FORTRAN_3D((*Ey)[mfi]),
- BL_TO_FORTRAN_3D((*Ez)[mfi]),
- BL_TO_FORTRAN_3D((*F)[mfi]),
- &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
- &WarpX::maxwell_fdtd_solver_id);
-
- }
-
- if (cost) {
- Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
- if (ipatch == 1) cbx.refine(rr);
- wt = (amrex::second() - wt) / cbx.d_numPts();
- (*cost)[mfi].plus(wt, cbx);
- }
- }
+ for ( MFIter mfi(*Ex,true); mfi.isValid(); ++mfi )
+ {
+ Real wt = amrex::second();
+
+ const Box& tex = mfi.tilebox(Ex_nodal_flag);
+ const Box& tey = mfi.tilebox(Ey_nodal_flag);
+ const Box& tez = mfi.tilebox(Ez_nodal_flag);
+
+ // Call picsar routine for each tile
+ warpx_push_evec(
+ tex.loVect(), tex.hiVect(),
+ tey.loVect(), tey.hiVect(),
+ tez.loVect(), tez.hiVect(),
+ BL_TO_FORTRAN_3D((*Ex)[mfi]),
+ BL_TO_FORTRAN_3D((*Ey)[mfi]),
+ BL_TO_FORTRAN_3D((*Ez)[mfi]),
+ BL_TO_FORTRAN_3D((*Bx)[mfi]),
+ BL_TO_FORTRAN_3D((*By)[mfi]),
+ BL_TO_FORTRAN_3D((*Bz)[mfi]),
+ BL_TO_FORTRAN_3D((*jx)[mfi]),
+ BL_TO_FORTRAN_3D((*jy)[mfi]),
+ BL_TO_FORTRAN_3D((*jz)[mfi]),
+ &mu_c2_dt,
+ &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
+
+ if (F) {
+ warpx_push_evec_f(
+ tex.loVect(), tex.hiVect(),
+ tey.loVect(), tey.hiVect(),
+ tez.loVect(), tez.hiVect(),
+ BL_TO_FORTRAN_3D((*Ex)[mfi]),
+ BL_TO_FORTRAN_3D((*Ey)[mfi]),
+ BL_TO_FORTRAN_3D((*Ez)[mfi]),
+ BL_TO_FORTRAN_3D((*F)[mfi]),
+ &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
+ &WarpX::maxwell_fdtd_solver_id);
+ }
+
+ if (cost) {
+ Box cbx = mfi.tilebox(IntVect{AMREX_D_DECL(0,0,0)});
+ if (patch_type == PatchType::coarse) cbx.refine(rr);
+ wt = (amrex::second() - wt) / cbx.d_numPts();
+ (*cost)[mfi].plus(wt, cbx);
+ }
}
- if (do_pml && pml[lev]->ok())
+ if (do_pml && pml[lev]->ok())
{
+ if (F) pml[lev]->ExchangeF(patch_type, F);
- for (int ipatch = 0; ipatch < npatches; ++ipatch)
- {
- const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
- const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
- const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
- int patch_level = (ipatch == 0) ? lev : lev-1;
- const std::array<Real,3>& dx = WarpX::CellSize(patch_level);
- const std::array<Real,3> dtsdx_c2 {c2dt/dx[0], c2dt/dx[1], c2dt/dx[2]};
+ const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
+ const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+ const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
#ifdef _OPENMP
#pragma omp parallel
#endif
- for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
+ for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
+ {
+ const Box& tex = mfi.tilebox(Ex_nodal_flag);
+ const Box& tey = mfi.tilebox(Ey_nodal_flag);
+ const Box& tez = mfi.tilebox(Ez_nodal_flag);
+
+ WRPX_PUSH_PML_EVEC(
+ tex.loVect(), tex.hiVect(),
+ tey.loVect(), tey.hiVect(),
+ tez.loVect(), tez.hiVect(),
+ BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
+ &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
+
+ if (pml_F)
{
- const Box& tex = mfi.tilebox(Ex_nodal_flag);
- const Box& tey = mfi.tilebox(Ey_nodal_flag);
- const Box& tez = mfi.tilebox(Ez_nodal_flag);
-
- WRPX_PUSH_PML_EVEC(
- tex.loVect(), tex.hiVect(),
- tey.loVect(), tey.hiVect(),
- tez.loVect(), tez.hiVect(),
- BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
- &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2]);
-
- if (pml_F)
- {
- WRPX_PUSH_PML_EVEC_F(
- tex.loVect(), tex.hiVect(),
- tey.loVect(), tey.hiVect(),
- tez.loVect(), tez.hiVect(),
- BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
- BL_TO_FORTRAN_3D((*pml_F )[mfi]),
- &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
- &WarpX::maxwell_fdtd_solver_id);
- }
+ WRPX_PUSH_PML_EVEC_F(
+ tex.loVect(), tex.hiVect(),
+ tey.loVect(), tey.hiVect(),
+ tez.loVect(), tez.hiVect(),
+ BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_F )[mfi]),
+ &dtsdx_c2[0], &dtsdx_c2[1], &dtsdx_c2[2],
+ &WarpX::maxwell_fdtd_solver_id);
}
}
}
}
+
void
WarpX::EvolveF (Real dt, DtType dt_type)
{
- if (!do_dive_cleaning) return;
+ if (!do_dive_cleaning) return;
- for (int lev = 0; lev <= finest_level; ++lev) {
- EvolveF(lev, dt, dt_type);
- }
+ for (int lev = 0; lev <= finest_level; ++lev) {
+ EvolveF(lev, dt, dt_type);
+ }
}
void
WarpX::EvolveF (int lev, Real dt, DtType dt_type)
{
- if (!do_dive_cleaning) return;
-
- BL_PROFILE("WarpX::EvolveF()");
+ if (!do_dive_cleaning) return;
- static constexpr Real c2inv = 1.0/(PhysConst::c*PhysConst::c);
- static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c;
-
- int npatches = (lev == 0) ? 1 : 2;
+ EvolveF(lev, PatchType::fine, dt, dt_type);
+ if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type);
+}
- for (int ipatch = 0; ipatch < npatches; ++ipatch)
- {
- int patch_level = (ipatch == 0) ? lev : lev-1;
- const auto& dx = WarpX::CellSize(patch_level);
- const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
+void
+WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type)
+{
+ if (!do_dive_cleaning) return;
- MultiFab *Ex, *Ey, *Ez, *rho, *F;
- if (ipatch == 0)
- {
- Ex = Efield_fp[lev][0].get();
- Ey = Efield_fp[lev][1].get();
- Ez = Efield_fp[lev][2].get();
- rho = rho_fp[lev].get();
- F = F_fp[lev].get();
- }
- else
- {
- Ex = Efield_cp[lev][0].get();
- Ey = Efield_cp[lev][1].get();
- Ez = Efield_cp[lev][2].get();
- rho = rho_cp[lev].get();
- F = F_cp[lev].get();
- }
+ BL_PROFILE("WarpX::EvolveF()");
- const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1;
+ static constexpr Real c2inv = 1.0/(PhysConst::c*PhysConst::c);
+ static constexpr Real mu_c2 = PhysConst::mu0*PhysConst::c*PhysConst::c;
- MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0);
- ComputeDivE(src, 0, {Ex,Ey,Ez}, dx);
- MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0);
- MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0);
+ int patch_level = (patch_type == PatchType::fine) ? lev : lev-1;
+ const auto& dx = WarpX::CellSize(patch_level);
+ const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]};
- if (do_pml && pml[lev]->ok())
- {
- const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
- const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+ MultiFab *Ex, *Ey, *Ez, *rho, *F;
+ if (patch_type == PatchType::fine)
+ {
+ Ex = Efield_fp[lev][0].get();
+ Ey = Efield_fp[lev][1].get();
+ Ez = Efield_fp[lev][2].get();
+ rho = rho_fp[lev].get();
+ F = F_fp[lev].get();
+ }
+ else
+ {
+ Ex = Efield_cp[lev][0].get();
+ Ey = Efield_cp[lev][1].get();
+ Ez = Efield_cp[lev][2].get();
+ rho = rho_cp[lev].get();
+ F = F_cp[lev].get();
+ }
+
+ const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1;
+
+ MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0);
+ ComputeDivE(src, 0, {Ex,Ey,Ez}, dx);
+ MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0);
+ MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0);
+
+ if (do_pml && pml[lev]->ok())
+ {
+ const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
+ const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
#ifdef _OPENMP
#pragma omp parallel
#endif
- for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi )
- {
- const Box& bx = mfi.tilebox();
- WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(),
- BL_TO_FORTRAN_ANYD((*pml_F )[mfi]),
- BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]),
- BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]),
- BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]),
- &dtsdx[0], &dtsdx[1], &dtsdx[2]);
- }
+ for ( MFIter mfi(*pml_F,true); mfi.isValid(); ++mfi )
+ {
+ const Box& bx = mfi.tilebox();
+ WRPX_PUSH_PML_F(bx.loVect(), bx.hiVect(),
+ BL_TO_FORTRAN_ANYD((*pml_F )[mfi]),
+ BL_TO_FORTRAN_ANYD((*pml_E[0])[mfi]),
+ BL_TO_FORTRAN_ANYD((*pml_E[1])[mfi]),
+ BL_TO_FORTRAN_ANYD((*pml_E[2])[mfi]),
+ &dtsdx[0], &dtsdx[1], &dtsdx[2]);
}
}
}
@@ -607,65 +601,65 @@ WarpX::EvolveF (int lev, Real dt, DtType dt_type)
void
WarpX::DampPML ()
{
- if (!do_pml) return;
-
- for (int lev = 0; lev <= finest_level; ++lev) {
- DampPML(lev);
- }
+ for (int lev = 0; lev <= finest_level; ++lev) {
+ DampPML(lev);
+ }
}
void
WarpX::DampPML (int lev)
{
- if (!do_pml) return;
+ DampPML(lev, PatchType::fine);
+ if (lev > 0) DampPML(lev, PatchType::coarse);
+}
- BL_PROFILE("WarpX::DampPML()");
+void
+WarpX::DampPML (int lev, PatchType patch_type)
+{
+ if (!do_pml) return;
- int npatches = (lev == 0) ? 1 : 2;
+ BL_PROFILE("WarpX::DampPML()");
- for (int ipatch = 0; ipatch < npatches; ++ipatch)
+ if (pml[lev]->ok())
{
- if (pml[lev]->ok())
- {
- const auto& pml_E = (ipatch==0) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
- const auto& pml_B = (ipatch==0) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
- const auto& pml_F = (ipatch==0) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
- const auto& sigba = (ipatch==0) ? pml[lev]->GetMultiSigmaBox_fp()
- : pml[lev]->GetMultiSigmaBox_cp();
+ const auto& pml_E = (patch_type == PatchType::fine) ? pml[lev]->GetE_fp() : pml[lev]->GetE_cp();
+ const auto& pml_B = (patch_type == PatchType::fine) ? pml[lev]->GetB_fp() : pml[lev]->GetB_cp();
+ const auto& pml_F = (patch_type == PatchType::fine) ? pml[lev]->GetF_fp() : pml[lev]->GetF_cp();
+ const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp()
+ : pml[lev]->GetMultiSigmaBox_cp();
#ifdef _OPENMP
#pragma omp parallel
#endif
- for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
- {
- const Box& tex = mfi.tilebox(Ex_nodal_flag);
- const Box& tey = mfi.tilebox(Ey_nodal_flag);
- const Box& tez = mfi.tilebox(Ez_nodal_flag);
- const Box& tbx = mfi.tilebox(Bx_nodal_flag);
- const Box& tby = mfi.tilebox(By_nodal_flag);
- const Box& tbz = mfi.tilebox(Bz_nodal_flag);
-
- WRPX_DAMP_PML(tex.loVect(), tex.hiVect(),
- tey.loVect(), tey.hiVect(),
- tez.loVect(), tez.hiVect(),
- tbx.loVect(), tbx.hiVect(),
- tby.loVect(), tby.hiVect(),
- tbz.loVect(), tbz.hiVect(),
- BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
- BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
- BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
- WRPX_PML_TO_FORTRAN(sigba[mfi]));
-
- if (pml_F) {
- const Box& tnd = mfi.nodaltilebox();
- WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(),
- BL_TO_FORTRAN_3D((*pml_F)[mfi]),
- WRPX_PML_TO_FORTRAN(sigba[mfi]));
- }
- }
+ for ( MFIter mfi(*pml_E[0],true); mfi.isValid(); ++mfi )
+ {
+ const Box& tex = mfi.tilebox(Ex_nodal_flag);
+ const Box& tey = mfi.tilebox(Ey_nodal_flag);
+ const Box& tez = mfi.tilebox(Ez_nodal_flag);
+ const Box& tbx = mfi.tilebox(Bx_nodal_flag);
+ const Box& tby = mfi.tilebox(By_nodal_flag);
+ const Box& tbz = mfi.tilebox(Bz_nodal_flag);
+
+ WRPX_DAMP_PML(tex.loVect(), tex.hiVect(),
+ tey.loVect(), tey.hiVect(),
+ tez.loVect(), tez.hiVect(),
+ tbx.loVect(), tbx.hiVect(),
+ tby.loVect(), tby.hiVect(),
+ tbz.loVect(), tbz.hiVect(),
+ BL_TO_FORTRAN_3D((*pml_E[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_E[2])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[0])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[1])[mfi]),
+ BL_TO_FORTRAN_3D((*pml_B[2])[mfi]),
+ WRPX_PML_TO_FORTRAN(sigba[mfi]));
+
+ if (pml_F) {
+ const Box& tnd = mfi.nodaltilebox();
+ WRPX_DAMP_PML_F(tnd.loVect(), tnd.hiVect(),
+ BL_TO_FORTRAN_3D((*pml_F)[mfi]),
+ WRPX_PML_TO_FORTRAN(sigba[mfi]));
+ }
}
}
}
diff --git a/Source/WarpXPML.H b/Source/WarpXPML.H
index c161ac596..d2ea26898 100644
--- a/Source/WarpXPML.H
+++ b/Source/WarpXPML.H
@@ -86,6 +86,8 @@ private:
amrex::Real dt_E = -1.e10;
};
+enum struct PatchType : int;
+
class PML
{
public:
@@ -113,13 +115,21 @@ public:
const std::array<amrex::MultiFab*,3>& B_cp);
void ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp,
const std::array<amrex::MultiFab*,3>& E_cp);
+ void ExchangeB (PatchType patch_type,
+ const std::array<amrex::MultiFab*,3>& Bp);
+ void ExchangeE (PatchType patch_type,
+ const std::array<amrex::MultiFab*,3>& Ep);
void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp);
+ void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp);
void FillBoundary ();
void FillBoundaryE ();
void FillBoundaryB ();
void FillBoundaryF ();
+ void FillBoundaryE (PatchType patch_type);
+ void FillBoundaryB (PatchType patch_type);
+ void FillBoundaryF (PatchType patch_type);
bool ok () const { return m_ok; }
diff --git a/Source/WarpXPML.cpp b/Source/WarpXPML.cpp
index cb8e8cc12..f9b924725 100644
--- a/Source/WarpXPML.cpp
+++ b/Source/WarpXPML.cpp
@@ -520,17 +520,25 @@ void
PML::ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp,
const std::array<amrex::MultiFab*,3>& B_cp)
{
- if (pml_B_fp[0])
+ ExchangeB(PatchType::fine, B_fp);
+ ExchangeB(PatchType::coarse, B_cp);
+}
+
+void
+PML::ExchangeB (PatchType patch_type,
+ const std::array<amrex::MultiFab*,3>& Bp)
+{
+ if (patch_type == PatchType::fine && pml_B_fp[0] && Bp[0])
{
- Exchange(*pml_B_fp[0], *B_fp[0], *m_geom);
- Exchange(*pml_B_fp[1], *B_fp[1], *m_geom);
- Exchange(*pml_B_fp[2], *B_fp[2], *m_geom);
+ Exchange(*pml_B_fp[0], *Bp[0], *m_geom);
+ Exchange(*pml_B_fp[1], *Bp[1], *m_geom);
+ Exchange(*pml_B_fp[2], *Bp[2], *m_geom);
}
- if (B_cp[0] && pml_B_cp[0])
+ else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0])
{
- Exchange(*pml_B_cp[0], *B_cp[0], *m_cgeom);
- Exchange(*pml_B_cp[1], *B_cp[1], *m_cgeom);
- Exchange(*pml_B_cp[2], *B_cp[2], *m_cgeom);
+ Exchange(*pml_B_cp[0], *Bp[0], *m_cgeom);
+ Exchange(*pml_B_cp[1], *Bp[1], *m_cgeom);
+ Exchange(*pml_B_cp[2], *Bp[2], *m_cgeom);
}
}
@@ -538,25 +546,43 @@ void
PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp,
const std::array<amrex::MultiFab*,3>& E_cp)
{
- if (pml_E_fp[0])
+ ExchangeB(PatchType::fine, E_fp);
+ ExchangeB(PatchType::coarse, E_cp);
+}
+
+void
+PML::ExchangeE (PatchType patch_type,
+ const std::array<amrex::MultiFab*,3>& Ep)
+{
+ if (patch_type == PatchType::fine && pml_E_fp[0] && Ep[0])
{
- Exchange(*pml_E_fp[0], *E_fp[0], *m_geom);
- Exchange(*pml_E_fp[1], *E_fp[1], *m_geom);
- Exchange(*pml_E_fp[2], *E_fp[2], *m_geom);
+ Exchange(*pml_E_fp[0], *Ep[0], *m_geom);
+ Exchange(*pml_E_fp[1], *Ep[1], *m_geom);
+ Exchange(*pml_E_fp[2], *Ep[2], *m_geom);
}
- if (E_cp[0] && pml_E_cp[0])
+ else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0])
{
- Exchange(*pml_E_cp[0], *E_cp[0], *m_cgeom);
- Exchange(*pml_E_cp[1], *E_cp[1], *m_cgeom);
- Exchange(*pml_E_cp[2], *E_cp[2], *m_cgeom);
+ Exchange(*pml_E_cp[0], *Ep[0], *m_cgeom);
+ Exchange(*pml_E_cp[1], *Ep[1], *m_cgeom);
+ Exchange(*pml_E_cp[2], *Ep[2], *m_cgeom);
}
}
void
PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp)
{
- if (pml_F_fp) Exchange(*pml_F_fp, *F_fp, *m_geom);
- if (pml_F_cp) Exchange(*pml_F_cp, *F_cp, *m_cgeom);
+ ExchangeF(PatchType::fine, F_fp);
+ ExchangeF(PatchType::coarse, F_cp);
+}
+
+void
+PML::ExchangeF (PatchType patch_type, MultiFab* Fp)
+{
+ if (patch_type == PatchType::fine && pml_F_fp && Fp) {
+ Exchange(*pml_F_fp, *Fp, *m_geom);
+ } else if (patch_type == PatchType::coarse && pml_F_cp && Fp) {
+ Exchange(*pml_F_cp, *Fp, *m_cgeom);
+ }
}
void
@@ -613,56 +639,74 @@ PML::FillBoundary ()
void
PML::FillBoundaryE ()
{
- if (pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0)
+ FillBoundaryE(PatchType::fine);
+ FillBoundaryE(PatchType::coarse);
+}
+
+void
+PML::FillBoundaryE (PatchType patch_type)
+{
+ if (patch_type == PatchType::fine && pml_E_fp[0] && pml_E_fp[0]->nGrowVect().max() > 0)
{
- const auto& period = m_geom->periodicity();
- pml_E_fp[0]->FillBoundary(period);
- pml_E_fp[1]->FillBoundary(period);
- pml_E_fp[2]->FillBoundary(period);
+ const auto& period = m_geom->periodicity();
+ pml_E_fp[0]->FillBoundary(period);
+ pml_E_fp[1]->FillBoundary(period);
+ pml_E_fp[2]->FillBoundary(period);
}
-
- if (pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0)
+ else if (patch_type == PatchType::coarse && pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0)
{
- const auto& period = m_cgeom->periodicity();
- pml_E_cp[0]->FillBoundary(period);
- pml_E_cp[1]->FillBoundary(period);
- pml_E_cp[2]->FillBoundary(period);
+ const auto& period = m_cgeom->periodicity();
+ pml_E_cp[0]->FillBoundary(period);
+ pml_E_cp[1]->FillBoundary(period);
+ pml_E_cp[2]->FillBoundary(period);
}
}
void
PML::FillBoundaryB ()
{
- if (pml_B_fp[0])
+ FillBoundaryB(PatchType::fine);
+ FillBoundaryB(PatchType::coarse);
+}
+
+void
+PML::FillBoundaryB (PatchType patch_type)
+{
+ if (patch_type == PatchType::fine && pml_B_fp[0])
{
- const auto& period = m_geom->periodicity();
- pml_B_fp[0]->FillBoundary(period);
- pml_B_fp[1]->FillBoundary(period);
- pml_B_fp[2]->FillBoundary(period);
+ const auto& period = m_geom->periodicity();
+ pml_B_fp[0]->FillBoundary(period);
+ pml_B_fp[1]->FillBoundary(period);
+ pml_B_fp[2]->FillBoundary(period);
}
-
- if (pml_B_cp[0])
+ else if (patch_type == PatchType::coarse && pml_B_cp[0])
{
- const auto& period = m_cgeom->periodicity();
- pml_B_cp[0]->FillBoundary(period);
- pml_B_cp[1]->FillBoundary(period);
- pml_B_cp[2]->FillBoundary(period);
+ const auto& period = m_cgeom->periodicity();
+ pml_B_cp[0]->FillBoundary(period);
+ pml_B_cp[1]->FillBoundary(period);
+ pml_B_cp[2]->FillBoundary(period);
}
}
void
PML::FillBoundaryF ()
{
- if (pml_F_fp && pml_F_fp->nGrowVect().max() > 0)
+ FillBoundaryF(PatchType::fine);
+ FillBoundaryF(PatchType::coarse);
+}
+
+void
+PML::FillBoundaryF (PatchType patch_type)
+{
+ if (patch_type == PatchType::fine && pml_F_fp && pml_F_fp->nGrowVect().max() > 0)
{
- const auto& period = m_geom->periodicity();
- pml_F_fp->FillBoundary(period);
+ const auto& period = m_geom->periodicity();
+ pml_F_fp->FillBoundary(period);
}
-
- if (pml_F_cp && pml_F_cp->nGrowVect().max() > 0)
+ else if (patch_type == PatchType::coarse && pml_F_cp && pml_F_cp->nGrowVect().max() > 0)
{
- const auto& period = m_cgeom->periodicity();
- pml_F_cp->FillBoundary(period);
+ const auto& period = m_cgeom->periodicity();
+ pml_F_cp->FillBoundary(period);
}
}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 12:12:56 +0000