aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Evolve/WarpXEvolveEM.cpp82
-rw-r--r--Source/Parallelization/GuardCellManager.H41
-rw-r--r--Source/Parallelization/GuardCellManager.cpp135
-rw-r--r--Source/Parallelization/Make.package2
-rw-r--r--Source/Parallelization/WarpXComm.cpp70
-rw-r--r--Source/Python/WarpXWrappers.cpp4
-rw-r--r--Source/WarpX.H37
-rw-r--r--Source/WarpX.cpp119
8 files changed, 317 insertions, 173 deletions
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp
index b5fd52bdc..8ad3b2401 100644
--- a/Source/Evolve/WarpXEvolveEM.cpp
+++ b/Source/Evolve/WarpXEvolveEM.cpp
@@ -75,8 +75,8 @@ WarpX::EvolveEM (int numsteps)
// Particles have p^{n} and x^{n}.
// is_synchronized is true.
if (is_synchronized) {
- FillBoundaryE();
- FillBoundaryB();
+ FillBoundaryE(guard_cells.ng_FieldGather);
+ FillBoundaryB(guard_cells.ng_FieldGather);
UpdateAuxilaryData();
// on first step, push p by -0.5*dt
for (int lev = 0; lev <= finest_level; ++lev) {
@@ -89,8 +89,10 @@ WarpX::EvolveEM (int numsteps)
} else {
// Beyond one step, we have E^{n} and B^{n}.
// Particles have p^{n-1/2} and x^{n}.
- FillBoundaryE();
- FillBoundaryB();
+ // This is probably overkill
+ FillBoundaryE(guard_cells.ngE);
+ // This is probably overkill
+ FillBoundaryB(guard_cells.ngE);
UpdateAuxilaryData();
}
@@ -185,8 +187,10 @@ WarpX::EvolveEM (int numsteps)
// slice gen //
if (to_make_plot || do_insitu || to_make_slice_plot)
{
- FillBoundaryE();
- FillBoundaryB();
+ // This is probably overkill
+ FillBoundaryE(guard_cells.ngE);
+ // This is probably overkill
+ FillBoundaryB(guard_cells.ngE);
UpdateAuxilaryData();
for (int lev = 0; lev <= finest_level; ++lev) {
@@ -237,8 +241,10 @@ WarpX::EvolveEM (int numsteps)
if (write_plot_file || do_insitu)
{
- FillBoundaryE();
- FillBoundaryB();
+ // This is probably overkill
+ FillBoundaryE(guard_cells.ngE);
+ // This is probably overkill
+ FillBoundaryB(guard_cells.ngE);
UpdateAuxilaryData();
for (int lev = 0; lev <= finest_level; ++lev) {
@@ -308,23 +314,21 @@ WarpX::OneStep_nosub (Real cur_time)
#ifdef WARPX_USE_PSATD
PushPSATD(dt[0]);
if (do_pml) DampPML();
- FillBoundaryE();
- FillBoundaryB();
+ FillBoundaryE(guard_cells.ngE);
+ FillBoundaryB(guard_cells.ngE);
#else
EvolveF(0.5*dt[0], DtType::FirstHalf);
- FillBoundaryF();
+ FillBoundaryF(guard_cells.ngF);
EvolveB(0.5*dt[0]); // We now have B^{n+1/2}
- FillBoundaryB();
+ FillBoundaryB(guard_cells.ngE_FieldSolver);
EvolveE(dt[0]); // We now have E^{n+1}
- FillBoundaryE();
+ FillBoundaryE(guard_cells.ngE_FieldSolver);
EvolveF(0.5*dt[0], DtType::SecondHalf);
EvolveB(0.5*dt[0]); // We now have B^{n+1}
if (do_pml) {
DampPML();
- FillBoundaryE();
}
- FillBoundaryB();
#endif
}
@@ -368,21 +372,21 @@ WarpX::OneStep_sub1 (Real curtime)
EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]);
EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf);
- FillBoundaryB(fine_lev, PatchType::fine);
- FillBoundaryF(fine_lev, PatchType::fine);
+ FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ngE);
+ FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ngF);
EvolveE(fine_lev, PatchType::fine, dt[fine_lev]);
- FillBoundaryE(fine_lev, PatchType::fine);
+ FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ngE);
EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]);
EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::SecondHalf);
if (do_pml) {
DampPML(fine_lev, PatchType::fine);
- FillBoundaryE(fine_lev, PatchType::fine);
+ FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ngE);
}
- FillBoundaryB(fine_lev, PatchType::fine);
+ FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ngE);
// ii) Push particles on the coarse patch and mother grid.
// Push the fields on the coarse patch and mother grid
@@ -394,19 +398,19 @@ WarpX::OneStep_sub1 (Real curtime)
EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]);
EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::FirstHalf);
- FillBoundaryB(fine_lev, PatchType::coarse);
- FillBoundaryF(fine_lev, PatchType::coarse);
+ FillBoundaryB(fine_lev, PatchType::coarse, guard_cells.ngE);
+ FillBoundaryF(fine_lev, PatchType::coarse, guard_cells.ngF);
EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]);
- FillBoundaryE(fine_lev, PatchType::coarse);
+ FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ngE);
EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]);
EvolveF(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev], DtType::FirstHalf);
- FillBoundaryB(coarse_lev, PatchType::fine);
- FillBoundaryF(coarse_lev, PatchType::fine);
+ FillBoundaryB(coarse_lev, PatchType::fine, guard_cells.ngE);
+ FillBoundaryF(coarse_lev, PatchType::fine, guard_cells.ngF);
EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]);
- FillBoundaryE(coarse_lev, PatchType::fine);
+ FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ngE);
// iii) Get auxiliary fields on the fine grid, at dt[fine_lev]
UpdateAuxilaryData();
@@ -422,22 +426,22 @@ WarpX::OneStep_sub1 (Real curtime)
EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]);
EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::FirstHalf);
- FillBoundaryB(fine_lev, PatchType::fine);
- FillBoundaryF(fine_lev, PatchType::fine);
+ FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ngE);
+ FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ngF);
EvolveE(fine_lev, PatchType::fine, dt[fine_lev]);
- FillBoundaryE(fine_lev, PatchType::fine);
+ FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ngE);
EvolveB(fine_lev, PatchType::fine, 0.5*dt[fine_lev]);
EvolveF(fine_lev, PatchType::fine, 0.5*dt[fine_lev], DtType::SecondHalf);
if (do_pml) {
DampPML(fine_lev, PatchType::fine);
- FillBoundaryE(fine_lev, PatchType::fine);
+ FillBoundaryE(fine_lev, PatchType::fine, guard_cells.ngE);
}
- FillBoundaryB(fine_lev, PatchType::fine);
- FillBoundaryF(fine_lev, PatchType::fine);
+ FillBoundaryB(fine_lev, PatchType::fine, guard_cells.ngE);
+ FillBoundaryF(fine_lev, PatchType::fine, guard_cells.ngF);
// v) Push the fields on the coarse patch and mother grid
// by only half a coarse step (second half)
@@ -446,7 +450,7 @@ WarpX::OneStep_sub1 (Real curtime)
AddRhoFromFineLevelandSumBoundary(coarse_lev, ncomps, ncomps);
EvolveE(fine_lev, PatchType::coarse, dt[fine_lev]);
- FillBoundaryE(fine_lev, PatchType::coarse);
+ FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ngE);
EvolveB(fine_lev, PatchType::coarse, dt[fine_lev]);
EvolveF(fine_lev, PatchType::coarse, dt[fine_lev], DtType::SecondHalf);
@@ -454,24 +458,24 @@ WarpX::OneStep_sub1 (Real curtime)
if (do_pml) {
DampPML(fine_lev, PatchType::coarse); // do it twice
DampPML(fine_lev, PatchType::coarse);
- FillBoundaryE(fine_lev, PatchType::coarse);
+ FillBoundaryE(fine_lev, PatchType::coarse, guard_cells.ngE);
}
- FillBoundaryB(fine_lev, PatchType::coarse);
- FillBoundaryF(fine_lev, PatchType::coarse);
+ FillBoundaryB(fine_lev, PatchType::coarse, guard_cells.ngE);
+ FillBoundaryF(fine_lev, PatchType::coarse, guard_cells.ngF);
EvolveE(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]);
- FillBoundaryE(coarse_lev, PatchType::fine);
+ FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ngE);
EvolveB(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev]);
EvolveF(coarse_lev, PatchType::fine, 0.5*dt[coarse_lev], DtType::SecondHalf);
if (do_pml) {
DampPML(coarse_lev, PatchType::fine);
- FillBoundaryE(coarse_lev, PatchType::fine);
+ FillBoundaryE(coarse_lev, PatchType::fine, guard_cells.ngE);
}
- FillBoundaryB(coarse_lev, PatchType::fine);
+ FillBoundaryB(coarse_lev, PatchType::fine, guard_cells.ngE);
}
void
diff --git a/Source/Parallelization/GuardCellManager.H b/Source/Parallelization/GuardCellManager.H
new file mode 100644
index 000000000..4b85a4332
--- /dev/null
+++ b/Source/Parallelization/GuardCellManager.H
@@ -0,0 +1,41 @@
+#ifndef GUARDCELLMANAGER_H_
+#define GUARDCELLMANAGER_H_
+
+#include <AMReX_IntVect.H>
+
+class guardCellManager{
+
+public:
+
+ void Init(
+ const bool do_subcycling,
+ const bool do_fdtd_nci_corr,
+ const bool do_nodal,
+ const bool do_moving_window,
+ const bool do_fft_mpi_dec,
+ const bool aux_is_nodal,
+ const int moving_window_dir,
+ const int nox,
+ const int nox_fft, const int noy_fft, const int noz_fft,
+ const int nci_corr_stencil,
+ const int maxwell_fdtd_solver_id,
+ const int max_level);
+
+ // Guard cells to initialize multifabs
+ amrex::IntVect ngExtra;
+ amrex::IntVect ngE;
+ amrex::IntVect ngJ;
+ amrex::IntVect ngRho;
+ amrex::IntVect ngF;
+ int ngF_int;
+
+ // Guard cells to exchange data
+ amrex::IntVect ngB_FieldSolver;
+ amrex::IntVect ngE_FieldSolver;
+ amrex::IntVect ng_FieldGather;
+ amrex::IntVect ngJ_CurrentDepo;
+ amrex::IntVect ng_MovingWindow;
+ amrex::IntVect ng_NCIFilter;
+};
+
+#endif // GUARDCELLMANAGER_H_
diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp
new file mode 100644
index 000000000..c790c3472
--- /dev/null
+++ b/Source/Parallelization/GuardCellManager.cpp
@@ -0,0 +1,135 @@
+#include "GuardCellManager.H"
+
+using namespace amrex;
+
+void
+guardCellManager::Init(
+ const bool do_subcycling,
+ const bool do_fdtd_nci_corr,
+ const bool do_nodal,
+ const bool do_moving_window,
+ const bool do_fft_mpi_dec,
+ const bool aux_is_nodal,
+ const int moving_window_dir,
+ const int nox,
+ const int nox_fft, const int noy_fft, const int noz_fft,
+ const int nci_corr_stencil,
+ const int maxwell_fdtd_solver_id,
+ const int max_level)
+{
+ // When using subcycling, the particles on the fine level perform two pushes
+ // before being redistributed ; therefore, we need one extra guard cell
+ // (the particles may move by 2*c*dt)
+ const int ngx_tmp = (max_level > 0 && do_subcycling == 1) ? nox+1 : nox;
+ const int ngy_tmp = (max_level > 0 && do_subcycling == 1) ? nox+1 : nox;
+ const int ngz_tmp = (max_level > 0 && do_subcycling == 1) ? nox+1 : nox;
+
+ // Ex, Ey, Ez, Bx, By, and Bz have the same number of ghost cells.
+ // jx, jy, jz and rho have the same number of ghost cells.
+ // E and B have the same number of ghost cells as j and rho if NCI filter is not used,
+ // but different number of ghost cells in z-direction if NCI filter is used.
+ // The number of cells should be even, in order to easily perform the
+ // interpolation from coarse grid to fine grid.
+ int ngx = (ngx_tmp % 2) ? ngx_tmp+1 : ngx_tmp; // Always even number
+ int ngy = (ngy_tmp % 2) ? ngy_tmp+1 : ngy_tmp; // Always even number
+ int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number
+ int ngz;
+ if (do_fdtd_nci_corr) {
+ int ng = ngz_tmp + nci_corr_stencil;
+ ngz = (ng % 2) ? ng+1 : ng;
+ } else {
+ ngz = ngz_nonci;
+ }
+
+ // J is only interpolated from fine to coarse (not coarse to fine)
+ // and therefore does not need to be even.
+ int ngJx = ngx_tmp;
+ int ngJy = ngy_tmp;
+ int ngJz = ngz_tmp;
+
+ // When calling the moving window (with one level of refinement), we shift
+ // the fine grid by 2 cells ; therefore, we need at least 2 guard cells
+ // on level 1. This may not be necessary for level 0.
+ if (do_moving_window) {
+ ngx = std::max(ngx,2);
+ ngy = std::max(ngy,2);
+ ngz = std::max(ngz,2);
+ ngJx = std::max(ngJx,2);
+ ngJy = std::max(ngJy,2);
+ ngJz = std::max(ngJz,2);
+ }
+
+/*
+#if (AMREX_SPACEDIM == 3)
+ IntVect ngE(ngx,ngy,ngz);
+ IntVect ngJ(ngJx,ngJy,ngJz);
+#elif (AMREX_SPACEDIM == 2)
+ IntVect ngE(ngx,ngz);
+ IntVect ngJ(ngJx,ngJz);
+#endif
+
+ IntVect ngRho = ngJ+1; //One extra ghost cell, so that it's safe to deposit charge density
+ // after pushing particle.
+ int ngF = (do_moving_window) ? 2 : 0;
+*/
+
+#if (AMREX_SPACEDIM == 3)
+ ngE = IntVect(ngx,ngy,ngz);
+ ngJ = IntVect(ngJx,ngJy,ngJz);
+#elif (AMREX_SPACEDIM == 2)
+ ngE = IntVect(ngx,ngz);
+ ngJ = IntVect(ngJx,ngJz);
+#endif
+
+ ngRho = ngJ+1; //One extra ghost cell, so that it's safe to deposit charge density
+ // after pushing particle.
+ ngF_int = (do_moving_window) ? 2 : 0;
+ // CKC solver requires one additional guard cell
+ if (maxwell_fdtd_solver_id == 1) ngF_int = std::max( ngF_int, 1 );
+ ngF = IntVect(AMREX_D_DECL(ngF_int, ngF_int, ngF_int));
+
+#ifdef WARPX_USE_PSATD
+ if (do_fft_mpi_dec == false){
+ // All boxes should have the same number of guard cells
+ // (to avoid temporary parallel copies)
+ // Thus take the max of the required number of guards for each field
+ // Also: the number of guard cell should be enough to contain
+ // the stencil of the FFT solver. Here, this number (`ngFFT`)
+ // is determined *empirically* to be the order of the solver
+ // for nodal, and half the order of the solver for staggered.
+ IntVect ngFFT;
+ if (do_nodal) {
+ ngFFT = IntVect(AMREX_D_DECL(nox_fft, noy_fft, noz_fft));
+ } else {
+ ngFFT = IntVect(AMREX_D_DECL(nox_fft/2, noy_fft/2, noz_fft/2));
+ }
+ for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){
+ int ng_required = ngFFT[i_dim];
+ // Get the max
+ ng_required = std::max( ng_required, ngE[i_dim] );
+ ng_required = std::max( ng_required, ngJ[i_dim] );
+ ng_required = std::max( ng_required, ngRho[i_dim] );
+ ng_required = std::max( ng_required, ngF[i_dim] );
+ // Set the guard cells to this max
+ ngE[i_dim] = ng_required;
+ ngJ[i_dim] = ng_required;
+ ngF[i_dim] = ng_required;
+ ngRho[i_dim] = ng_required;
+ ngF_int = ng_required;
+ }
+ }
+ ngF = IntVect(AMREX_D_DECL(ngF_int, ngF_int, ngF_int));
+#endif
+
+ ngExtra = IntVect(static_cast<int>(aux_is_nodal and !do_nodal));
+
+ // Guard cells for field solver
+ ngB_FieldSolver = IntVect(AMREX_D_DECL(1,1,1));
+ ngE_FieldSolver = IntVect(AMREX_D_DECL(1,1,1));
+ ng_MovingWindow = IntVect(AMREX_D_DECL(0,0,0)); // Multiplied by number of cells moved at each timestep
+ ng_MovingWindow[moving_window_dir] = 1;
+ int FGcell[4] = {0,1,1,2}; // Index is nox
+ ng_FieldGather = IntVect(AMREX_D_DECL(FGcell[nox],FGcell[nox],FGcell[nox]));
+ ngJ_CurrentDepo = ng_FieldGather;
+ ng_NCIFilter = IntVect(AMREX_D_DECL(0,0,4));
+}
diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package
index c74583522..cabce2da8 100644
--- a/Source/Parallelization/Make.package
+++ b/Source/Parallelization/Make.package
@@ -1,7 +1,9 @@
CEXE_sources += WarpXComm.cpp
CEXE_sources += WarpXRegrid.cpp
+CEXE_sources += GuardCellManager.cpp
CEXE_headers += WarpXSumGuardCells.H
CEXE_headers += WarpXComm_K.H
+CEXE_headers += GuardCellManager.H
INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parallelization
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 52df3dc25..0dae38e2e 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -321,41 +321,41 @@ WarpX::UpdateAuxilaryDataSameType ()
}
void
-WarpX::FillBoundaryB ()
+WarpX::FillBoundaryB (IntVect ng)
{
for (int lev = 0; lev <= finest_level; ++lev)
{
- FillBoundaryB(lev);
+ FillBoundaryB(lev, ng);
}
}
void
-WarpX::FillBoundaryE ()
+WarpX::FillBoundaryE (IntVect ng)
{
for (int lev = 0; lev <= finest_level; ++lev)
{
- FillBoundaryE(lev);
+ FillBoundaryE(lev, ng);
}
}
void
-WarpX::FillBoundaryF ()
+WarpX::FillBoundaryF (IntVect ng)
{
for (int lev = 0; lev <= finest_level; ++lev)
{
- FillBoundaryF(lev);
+ FillBoundaryF(lev, ng);
}
}
void
-WarpX::FillBoundaryE(int lev)
+WarpX::FillBoundaryE(int lev, IntVect ng)
{
- FillBoundaryE(lev, PatchType::fine);
- if (lev > 0) FillBoundaryE(lev, PatchType::coarse);
+ FillBoundaryE(lev, PatchType::fine, ng);
+ if (lev > 0) FillBoundaryE(lev, PatchType::coarse, ng);
}
void
-WarpX::FillBoundaryE (int lev, PatchType patch_type)
+WarpX::FillBoundaryE (int lev, PatchType patch_type, IntVect ng)
{
if (patch_type == PatchType::fine)
{
@@ -370,8 +370,11 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type)
}
const auto& period = Geom(lev).periodicity();
- Vector<MultiFab*> mf{Efield_fp[lev][0].get(),Efield_fp[lev][1].get(),Efield_fp[lev][2].get()};
- amrex::FillBoundary(mf, period);
+ Efield_fp[lev][0]->FillBoundary(0, Efield_fp[lev][0]->nComp(), ng, period);
+ Efield_fp[lev][1]->FillBoundary(0, Efield_fp[lev][1]->nComp(), ng, period);
+ Efield_fp[lev][2]->FillBoundary(0, Efield_fp[lev][2]->nComp(), ng, period);
+ // Vector<MultiFab*> mf{Efield_fp[lev][0].get(),Efield_fp[lev][1].get(),Efield_fp[lev][2].get()};
+ // amrex::FillBoundary(mf, period);
}
else if (patch_type == PatchType::coarse)
{
@@ -386,20 +389,23 @@ WarpX::FillBoundaryE (int lev, PatchType patch_type)
}
const auto& cperiod = Geom(lev-1).periodicity();
- Vector<MultiFab*> mf{Efield_cp[lev][0].get(),Efield_cp[lev][1].get(),Efield_cp[lev][2].get()};
- amrex::FillBoundary(mf, cperiod);
+ Efield_cp[lev][0]->FillBoundary(0, Efield_cp[lev][0]->nComp(), ng, cperiod);
+ Efield_cp[lev][1]->FillBoundary(0, Efield_cp[lev][1]->nComp(), ng, cperiod);
+ Efield_cp[lev][2]->FillBoundary(0, Efield_cp[lev][2]->nComp(), ng, cperiod);
+ // Vector<MultiFab*> mf{Efield_cp[lev][0].get(),Efield_cp[lev][1].get(),Efield_cp[lev][2].get()};
+ // amrex::FillBoundary(mf, cperiod);
}
}
void
-WarpX::FillBoundaryB (int lev)
+WarpX::FillBoundaryB (int lev, IntVect ng)
{
- FillBoundaryB(lev, PatchType::fine);
- if (lev > 0) FillBoundaryB(lev, PatchType::coarse);
+ FillBoundaryB(lev, PatchType::fine, ng);
+ if (lev > 0) FillBoundaryB(lev, PatchType::coarse, ng);
}
void
-WarpX::FillBoundaryB (int lev, PatchType patch_type)
+WarpX::FillBoundaryB (int lev, PatchType patch_type, IntVect ng)
{
if (patch_type == PatchType::fine)
{
@@ -413,8 +419,11 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type)
pml[lev]->FillBoundaryB(patch_type);
}
const auto& period = Geom(lev).periodicity();
- Vector<MultiFab*> mf{Bfield_fp[lev][0].get(),Bfield_fp[lev][1].get(),Bfield_fp[lev][2].get()};
- amrex::FillBoundary(mf, period);
+ Bfield_fp[lev][0]->FillBoundary(0, Bfield_fp[lev][0]->nComp(), ng, period);
+ Bfield_fp[lev][1]->FillBoundary(0, Bfield_fp[lev][1]->nComp(), ng, period);
+ Bfield_fp[lev][2]->FillBoundary(0, Bfield_fp[lev][2]->nComp(), ng, period);
+ // Vector<MultiFab*> mf{Bfield_fp[lev][0].get(),Bfield_fp[lev][1].get(),Bfield_fp[lev][2].get()};
+ // amrex::FillBoundary(mf, period);
}
else if (patch_type == PatchType::coarse)
{
@@ -428,20 +437,23 @@ WarpX::FillBoundaryB (int lev, PatchType patch_type)
pml[lev]->FillBoundaryB(patch_type);
}
const auto& cperiod = Geom(lev-1).periodicity();
- Vector<MultiFab*> mf{Bfield_cp[lev][0].get(),Bfield_cp[lev][1].get(),Bfield_cp[lev][2].get()};
- amrex::FillBoundary(mf, cperiod);
+ Bfield_cp[lev][0]->FillBoundary(0, Bfield_cp[lev][0]->nComp(), ng, cperiod);
+ Bfield_cp[lev][1]->FillBoundary(0, Bfield_cp[lev][1]->nComp(), ng, cperiod);
+ Bfield_cp[lev][2]->FillBoundary(0, Bfield_cp[lev][2]->nComp(), ng, cperiod);
+ // Vector<MultiFab*> mf{Bfield_cp[lev][0].get(),Bfield_cp[lev][1].get(),Bfield_cp[lev][2].get()};
+ // amrex::FillBoundary(mf, cperiod);
}
}
void
-WarpX::FillBoundaryF (int lev)
+WarpX::FillBoundaryF (int lev, IntVect ng)
{
- FillBoundaryF(lev, PatchType::fine);
- if (lev > 0) FillBoundaryF(lev, PatchType::coarse);
+ FillBoundaryF(lev, PatchType::fine, ng);
+ if (lev > 0) FillBoundaryF(lev, PatchType::coarse, ng);
}
void
-WarpX::FillBoundaryF (int lev, PatchType patch_type)
+WarpX::FillBoundaryF (int lev, PatchType patch_type, IntVect ng)
{
if (patch_type == PatchType::fine && F_fp[lev])
{
@@ -453,7 +465,8 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type)
}
const auto& period = Geom(lev).periodicity();
- F_fp[lev]->FillBoundary(period);
+ F_fp[lev]->FillBoundary(0, F_fp[lev]->nComp(), ng, period);
+ // F_fp[lev]->FillBoundary(period);
}
else if (patch_type == PatchType::coarse && F_cp[lev])
{
@@ -465,7 +478,8 @@ WarpX::FillBoundaryF (int lev, PatchType patch_type)
}
const auto& cperiod = Geom(lev-1).periodicity();
- F_cp[lev]->FillBoundary(cperiod);
+ F_cp[lev]->FillBoundary(0, F_cp[lev]->nComp(), ng, cperiod);
+ // F_cp[lev]->FillBoundary(cperiod);
}
}
diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp
index be9ec9519..3635bcd45 100644
--- a/Source/Python/WarpXWrappers.cpp
+++ b/Source/Python/WarpXWrappers.cpp
@@ -378,11 +378,11 @@ extern "C"
}
void warpx_FillBoundaryE () {
WarpX& warpx = WarpX::GetInstance();
- warpx.FillBoundaryE ();
+ warpx.FillBoundaryE (warpx.getngE());
}
void warpx_FillBoundaryB () {
WarpX& warpx = WarpX::GetInstance();
- warpx.FillBoundaryB ();
+ warpx.FillBoundaryB (warpx.getngE());
}
void warpx_SyncCurrent () {
WarpX& warpx = WarpX::GetInstance();
diff --git a/Source/WarpX.H b/Source/WarpX.H
index f55670cfb..da9431d32 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -25,6 +25,8 @@
#include <AMReX_Interpolater.H>
#include <AMReX_FillPatchUtil.H>
+#include "GuardCellManager.H"
+
#ifdef _OPENMP
# include <omp.h>
#endif
@@ -213,12 +215,12 @@ public:
void UpdateAuxilaryDataSameType ();
// Fill boundary cells including coarse/fine boundaries
- void FillBoundaryB ();
- void FillBoundaryE ();
- void FillBoundaryF ();
- void FillBoundaryE (int lev);
- void FillBoundaryB (int lev);
- void FillBoundaryF (int lev);
+ void FillBoundaryB (amrex::IntVect ng);
+ void FillBoundaryE (amrex::IntVect ng);
+ void FillBoundaryF (amrex::IntVect ng);
+ void FillBoundaryE (int lev, amrex::IntVect ng);
+ void FillBoundaryB (int lev, amrex::IntVect ng);
+ void FillBoundaryF (int lev, amrex::IntVect ng);
void SyncCurrent ();
void SyncRho ();
@@ -298,6 +300,9 @@ public:
const std::array<const amrex::MultiFab*, 3>& B,
const std::array<amrex::Real,3>& dx, int ngrow);
+ const amrex::IntVect getngE() const { return guard_cells.ngE; };
+ const amrex::IntVect getngF() const { return guard_cells.ngF; };
+
protected:
//! Tagging cells for refinement
@@ -335,9 +340,9 @@ private:
///
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 FillBoundaryB (int lev, PatchType patch_type, amrex::IntVect ng);
+ void FillBoundaryE (int lev, PatchType patch_type, amrex::IntVect ng);
+ void FillBoundaryF (int lev, PatchType patch_type, amrex::IntVect ng);
void OneStep_nosub (amrex::Real t);
void OneStep_sub1 (amrex::Real t);
@@ -470,7 +475,8 @@ private:
void AllocLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
const amrex::IntVect& ngE, const amrex::IntVect& ngJ,
- const amrex::IntVect& ngRho, int ngF);
+ const amrex::IntVect& ngRho, int ngF, const amrex::IntVect& ngextra,
+ const bool aux_is_nodal);
amrex::Vector<int> istep; // which step?
amrex::Vector<int> nsubsteps; // how many substeps on each level?
@@ -603,6 +609,8 @@ private:
bool is_synchronized = true;
+ guardCellManager guard_cells;
+
//Slice Parameters
int slice_max_grid_size;
int slice_plot_int = -1;
@@ -628,18 +636,19 @@ private:
amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_cp_fft;
#endif
+ bool fft_hybrid_mpi_decomposition = false;
+ int nox_fft = 16;
+ int noy_fft = 16;
+ int noz_fft = 16;
+
#ifdef WARPX_USE_PSATD
private:
void EvolvePSATD (int numsteps);
void PushPSATD (amrex::Real dt);
void PushPSATD_localFFT (int lev, amrex::Real dt);
- bool fft_hybrid_mpi_decomposition = false;
int ngroups_fft = 4;
int fftw_plan_measure = 1;
- int nox_fft = 16;
- int noy_fft = 16;
- int noz_fft = 16;
amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_fp;
amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_cp;
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index eef033236..1915434b6 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -700,58 +700,34 @@ WarpX::ClearLevel (int lev)
void
WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& dm)
{
- // When using subcycling, the particles on the fine level perform two pushes
- // before being redistributed ; therefore, we need one extra guard cell
- // (the particles may move by 2*c*dt)
- const int ngx_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::nox+1 : WarpX::nox;
- const int ngy_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noy+1 : WarpX::noy;
- const int ngz_tmp = (maxLevel() > 0 && do_subcycling == 1) ? WarpX::noz+1 : WarpX::noz;
-
- // Ex, Ey, Ez, Bx, By, and Bz have the same number of ghost cells.
- // jx, jy, jz and rho have the same number of ghost cells.
- // E and B have the same number of ghost cells as j and rho if NCI filter is not used,
- // but different number of ghost cells in z-direction if NCI filter is used.
- // The number of cells should be even, in order to easily perform the
- // interpolation from coarse grid to fine grid.
- int ngx = (ngx_tmp % 2) ? ngx_tmp+1 : ngx_tmp; // Always even number
- int ngy = (ngy_tmp % 2) ? ngy_tmp+1 : ngy_tmp; // Always even number
- int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number
- int ngz;
- if (WarpX::use_fdtd_nci_corr) {
- int ng = ngz_tmp + NCIGodfreyFilter::m_stencil_width;
- ngz = (ng % 2) ? ng+1 : ng;
- } else {
- ngz = ngz_nonci;
- }
- // J is only interpolated from fine to coarse (not coarse to fine)
- // and therefore does not need to be even.
- int ngJx = ngx_tmp;
- int ngJy = ngy_tmp;
- int ngJz = ngz_tmp;
-
- // When calling the moving window (with one level of refinement), we shift
- // the fine grid by 2 cells ; therefore, we need at least 2 guard cells
- // on level 1. This may not be necessary for level 0.
- if (do_moving_window) {
- ngx = std::max(ngx,2);
- ngy = std::max(ngy,2);
- ngz = std::max(ngz,2);
- ngJx = std::max(ngJx,2);
- ngJy = std::max(ngJy,2);
- ngJz = std::max(ngJz,2);
- }
-
-#if (AMREX_SPACEDIM == 3)
- IntVect ngE(ngx,ngy,ngz);
- IntVect ngJ(ngJx,ngJy,ngJz);
-#elif (AMREX_SPACEDIM == 2)
- IntVect ngE(ngx,ngz);
- IntVect ngJ(ngJx,ngJz);
-#endif
+ bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving);
- IntVect ngRho = ngJ+1; //One extra ghost cell, so that it's safe to deposit charge density
- // after pushing particle.
+ guard_cells.Init(
+ do_subcycling,
+ WarpX::use_fdtd_nci_corr,
+ do_nodal,
+ do_moving_window,
+ fft_hybrid_mpi_decomposition,
+ aux_is_nodal,
+ moving_window_dir,
+ WarpX::nox,
+ nox_fft, noy_fft, noz_fft,
+ NCIGodfreyFilter::m_stencil_width,
+ maxwell_fdtd_solver_id,
+ maxLevel());
+
+ IntVect ngE = guard_cells.ngE;
+ IntVect ngJ = guard_cells.ngJ;
+ IntVect ngRho = guard_cells.ngRho;
+ int ngF_int = guard_cells.ngF_int;
+ IntVect ngextra = guard_cells.ngExtra;
+
+ Print()<<"ngE "<<ngE<<'\n';
+ Print()<<"ngJ "<<ngJ<<'\n';
+ Print()<<"ngRho "<<ngRho<<'\n';
+ Print()<<"ngF_int "<<ngF_int<<'\n';
+ Print()<<"ngextra "<<ngextra<<'\n';
if (mypc->nSpeciesDepositOnMainGrid() && n_current_deposition_buffer == 0) {
n_current_deposition_buffer = 1;
@@ -769,47 +745,13 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d
n_field_gather_buffer = n_current_deposition_buffer + 1;
}
- int ngF = (do_moving_window) ? 2 : 0;
- // CKC solver requires one additional guard cell
- if (maxwell_fdtd_solver_id == 1) ngF = std::max( ngF, 1 );
-
-#ifdef WARPX_USE_PSATD
- if (fft_hybrid_mpi_decomposition == false){
- // All boxes should have the same number of guard cells
- // (to avoid temporary parallel copies)
- // Thus take the max of the required number of guards for each field
- // Also: the number of guard cell should be enough to contain
- // the stencil of the FFT solver. Here, this number (`ngFFT`)
- // is determined *empirically* to be the order of the solver
- // for nodal, and half the order of the solver for staggered.
- IntVect ngFFT;
- if (do_nodal) {
- ngFFT = IntVect(AMREX_D_DECL(nox_fft, noy_fft, noz_fft));
- } else {
- ngFFT = IntVect(AMREX_D_DECL(nox_fft/2, noy_fft/2, noz_fft/2));
- }
- for (int i_dim=0; i_dim<AMREX_SPACEDIM; i_dim++ ){
- int ng_required = ngFFT[i_dim];
- // Get the max
- ng_required = std::max( ng_required, ngE[i_dim] );
- ng_required = std::max( ng_required, ngJ[i_dim] );
- ng_required = std::max( ng_required, ngRho[i_dim] );
- ng_required = std::max( ng_required, ngF );
- // Set the guard cells to this max
- ngE[i_dim] = ng_required;
- ngJ[i_dim] = ng_required;
- ngRho[i_dim] = ng_required;
- ngF = ng_required;
- }
- }
-#endif
-
- AllocLevelMFs(lev, ba, dm, ngE, ngJ, ngRho, ngF);
+ AllocLevelMFs(lev, ba, dm, ngE, ngJ, ngRho, ngF_int, ngextra, aux_is_nodal);
}
void
WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm,
- const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, int ngF)
+ const IntVect& ngE, const IntVect& ngJ, const IntVect& ngRho, int ngF,
+ const IntVect& ngextra, const bool aux_is_nodal)
{
#if defined WARPX_DIM_RZ
@@ -821,9 +763,6 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
ncomps = n_rz_azimuthal_modes*2 - 1;
#endif
- bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving);
- IntVect ngextra(static_cast<int>(aux_is_nodal and !do_nodal));
-
//
// The fine patch
//
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-14 01:43:22 +0000