From ab7f93834aa20999c9abb402df65309748339551 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 3 May 2019 11:29:31 -0700 Subject: Increase number of guard cells for the PML --- Source/BoundaryConditions/PML.cpp | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index c3449cecd..d8052c4fa 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -343,10 +343,30 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, DistributionMapping dm{ba}; - int nge = 2; - int ngb = 2; - int ngf = (do_moving_window) ? 2 : 0; - if (WarpX::maxwell_fdtd_solver_id == 1) ngf = std::max( ngf, 1 ); + // Define the number of guard cells in each direction, for E, B, and F + IntVect nge = IntVect(AMREX_D_DECL(2, 2, 2)); + IntVect ngb = IntVect(AMREX_D_DECL(2, 2, 2)); + int ngf_int = (do_moving_window) ? 2 : 0; + if (WarpX::maxwell_fdtd_solver_id == 1) ngf_int = std::max( ngf_int, 1 ); + IntVect ngf = IntVect(AMREX_D_DECL(ngf_int, ngf_int, ngf_int)); +#ifdef WARPX_USE_PSATD + // Increase the number of guard cells, in order to fit the extent + // of the stencil for the spectral solver + 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)); + } + // Set the number of guard cells to the maximum of each field + // (all fields should have the same number of guard cells) + ngFFT = ngFFT.max(nge); + ngFFT = ngFFT.max(ngb); + ngFFT = ngFFT.max(ngf); + nge = ngFFT; + ngb = ngFFT; + ngf = ngFFT; + #endif pml_E_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::Ex_nodal_flag), dm, 3, nge)); pml_E_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::Ey_nodal_flag), dm, 3, nge)); @@ -372,9 +392,10 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, if (cgeom) { - - nge = 1; - ngb = 1; +#ifndef WARPX_USE_PSATD + nge = IntVect(AMREX_D_DECL(1, 1, 1)); + ngb = IntVect(AMREX_D_DECL(1, 1, 1)); +#endif BoxArray grid_cba = grid_ba; grid_cba.coarsen(ref_ratio); -- cgit v1.2.3 From 257e71c2eaca05416122e8b748c13a324e82d831 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Fri, 3 May 2019 15:45:09 -0700 Subject: Prepare spectralsolver for PML --- Source/BoundaryConditions/PML.H | 9 +++++++++ Source/BoundaryConditions/PML.cpp | 1 - Source/FieldSolver/SpectralSolver/SpectralSolver.H | 3 ++- Source/FieldSolver/SpectralSolver/SpectralSolver.cpp | 20 ++++++++++++++++---- 4 files changed, 27 insertions(+), 6 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 25dfc7996..39dfa5a25 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -6,6 +6,10 @@ #include #include +#ifdef WARPX_USE_PSATD +#include +#endif + #if (AMREX_SPACEDIM == 3) #define WRPX_PML_TO_FORTRAN(x) \ @@ -161,6 +165,11 @@ private: std::unique_ptr sigba_fp; std::unique_ptr sigba_cp; +#ifdef WARPX_USE_PSATD + std::unique_ptr spectral_solver_fp; + std::unique_ptr spectral_solver_cp; +#endif + static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell); diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index d8052c4fa..cac43faa5 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -425,7 +425,6 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); } - } BoxArray diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index d4019a9a3..c570b017b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -23,7 +23,8 @@ class SpectralSolver const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ); + const amrex::RealVect dx, const amrex::Real dt, + const bool pml=false ); /* \brief Transform the component `i_comp` of MultiFab `mf` * to spectral space, and store the corresponding result internally diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index a91fcbc47..80555a7b3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -7,13 +7,22 @@ * This function selects the spectral algorithm to be used, allocates the * corresponding coefficients for the discretized field update equation, * and prepares the structures that store the fields in spectral space. + * + * \param norder_x Order of accuracy of the spatial derivatives along x + * \param norder_y Order of accuracy of the spatial derivatives along y + * \param norder_z Order of accuracy of the spatial derivatives along z + * \param nodal Whether the solver is applied to a nodal or staggered grid + * \param dx Cell size along each dimension + * \param dt Time step + * \param pml Whether the boxes in which the solver is applied are PML boxes */ SpectralSolver::SpectralSolver( const amrex::BoxArray& realspace_ba, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, const int norder_z, const bool nodal, - const amrex::RealVect dx, const amrex::Real dt ) { + const amrex::RealVect dx, const amrex::Real dt, + const bool pml ) { // Initialize all structures using the same distribution mapping dm @@ -24,10 +33,13 @@ SpectralSolver::SpectralSolver( // - Select the algorithm depending on the input parameters // Initialize the corresponding coefficients over k space - // TODO: Add more algorithms + selection depending on input parameters - // For the moment, this only uses the standard PsatdAlgorithm - algorithm = std::unique_ptr( new PsatdAlgorithm( + if (pml) { + algorithm = std::unique_ptr( new PMLPsatdAlgorithm( + k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + } else { + algorithm = std::unique_ptr( new PsatdAlgorithm( k_space, dm, norder_x, norder_y, norder_z, nodal, dt ) ); + } // - Initialize arrays for fields in spectral space + FFT plans field_data = SpectralFieldData( realspace_ba, k_space, dm, -- cgit v1.2.3 From f71b9c0c409ebf1750b9c40e861bfeefac4247ca Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 4 May 2019 17:57:14 -0700 Subject: Initialize spectral solver --- Source/BoundaryConditions/PML.H | 4 +++- Source/BoundaryConditions/PML.cpp | 18 +++++++++++++++++- .../FieldSolver/SpectralSolver/SpectralSolver.cpp | 1 + Source/Initialization/WarpXInitData.cpp | 21 ++++++++++++--------- 4 files changed, 33 insertions(+), 11 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 39dfa5a25..15c70c10b 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -104,7 +104,9 @@ class PML public: PML (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const amrex::Geometry* geom, const amrex::Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window); + int ncell, int delta, int ref_ratio, amrex::Real dt, + int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int do_dive_cleaning, int do_moving_window); void ComputePMLFactors (amrex::Real dt); diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index cac43faa5..831339a2f 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -329,7 +329,9 @@ MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, const Geometry* geom, const Geometry* cgeom, - int ncell, int delta, int ref_ratio, int do_dive_cleaning, int do_moving_window) + int ncell, int delta, int ref_ratio, Real dt, + int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int do_dive_cleaning, int do_moving_window) : m_geom(geom), m_cgeom(cgeom) { @@ -390,6 +392,13 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); +#ifdef WARPX_USE_PSATD + const bool in_pml = true; // Tells spectral solver to use split-PML equations + const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))}; + spectral_solver_fp.reset( new SpectralSolver( ba, dm, + nox_fft, noy_fft, noz_fft, do_nodal, dx, dt, in_pml ) ); +#endif + if (cgeom) { #ifndef WARPX_USE_PSATD @@ -424,6 +433,13 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, } sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + + #ifdef WARPX_USE_PSATD + const bool in_pml = true; // Tells spectral solver to use split-PML equations + const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))}; + spectral_solver_fp.reset( new SpectralSolver( cba, cdm, + nox_fft, noy_fft, noz_fft, do_nodal, cdx, dt, in_pml ) ); + #endif } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index 80555a7b3..4b9def013 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -1,6 +1,7 @@ #include #include #include +#include /* \brief Initialize the spectral Maxwell solver * diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 23637ec97..f59f5d81f 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -87,7 +87,7 @@ WarpX::InitDiagnostics () { const Real* current_lo = geom[0].ProbLo(); const Real* current_hi = geom[0].ProbHi(); Real dt_boost = dt[0]; - + // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0 Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); @@ -96,7 +96,7 @@ WarpX::InitDiagnostics () { zmax_lab, moving_window_v, dt_snapshots_lab, num_snapshots_lab, gamma_boost, - t_new[0], dt_boost, + t_new[0], dt_boost, moving_window_dir, geom[0])); } } @@ -117,10 +117,10 @@ WarpX::InitFromScratch () InitPML(); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC if (do_electrostatic) { getLevelMasks(masks); - + // the plus one is to convert from num_cells to num_nodes getLevelMasks(gather_masks, n_buffer + 1); } @@ -133,13 +133,16 @@ WarpX::InitPML () if (do_pml) { pml[0].reset(new PML(boxArray(0), DistributionMap(0), &Geom(0), nullptr, - pml_ncell, pml_delta, 0, do_dive_cleaning, do_moving_window)); + pml_ncell, pml_delta, 0, dt[0], + nox_fft, noy_fft, noz_fft, do_nodal, + do_dive_cleaning, do_moving_window)); for (int lev = 1; lev <= finest_level; ++lev) { pml[lev].reset(new PML(boxArray(lev), DistributionMap(lev), &Geom(lev), &Geom(lev-1), - pml_ncell, pml_delta, refRatio(lev-1)[0], do_dive_cleaning, - do_moving_window)); + pml_ncell, pml_delta, refRatio(lev-1)[0], dt[lev], + nox_fft, noy_fft, noz_fft, do_nodal, + do_dive_cleaning, do_moving_window)); } } } @@ -222,7 +225,7 @@ WarpX::InitOpenbc () Vector alllohi(6*nprocs,100000); MPI_Allgather(lohi, 6, MPI_INT, alllohi.data(), 6, MPI_INT, ParallelDescriptor::Communicator()); - + BoxList bl{IndexType::TheNodeType()}; for (int i = 0; i < nprocs; ++i) { @@ -248,7 +251,7 @@ WarpX::InitOpenbc () rho_openbc.copy(*rho, 0, 0, 1, rho->nGrow(), 0, gm.periodicity(), FabArrayBase::ADD); const Real* dx = gm.CellSize(); - + warpx_openbc_potential(rho_openbc[myproc].dataPtr(), phi_openbc[myproc].dataPtr(), dx); BoxArray nba = boxArray(lev); -- cgit v1.2.3 From 39849af4464a48f2c7420d2ae94b1f9a261c9309 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sat, 4 May 2019 21:46:11 -0700 Subject: Perform transform --- Source/BoundaryConditions/PML.H | 4 ++++ Source/BoundaryConditions/PML.cpp | 50 ++++++++++++++++++++++++++++++++++----- Source/FieldSolver/WarpXFFT.cpp | 9 ++++++- 3 files changed, 56 insertions(+), 7 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 15c70c10b..ee38a5af9 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -124,6 +124,10 @@ public: const MultiSigmaBox& GetMultiSigmaBox_cp () const { return *sigba_cp; } +#ifdef WARPX_USE_PSATD + void PushPSATD (); +#endif + void ExchangeB (const std::array& B_fp, const std::array& B_cp); void ExchangeE (const std::array& E_fp, diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 831339a2f..77693915e 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -434,12 +434,12 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); - #ifdef WARPX_USE_PSATD - const bool in_pml = true; // Tells spectral solver to use split-PML equations - const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))}; - spectral_solver_fp.reset( new SpectralSolver( cba, cdm, - nox_fft, noy_fft, noz_fft, do_nodal, cdx, dt, in_pml ) ); - #endif +#ifdef WARPX_USE_PSATD + const bool in_pml = true; // Tells spectral solver to use split-PML equations + const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))}; + spectral_solver_fp.reset( new SpectralSolver( cba, cdm, + nox_fft, noy_fft, noz_fft, do_nodal, cdx, dt, in_pml ) ); +#endif } } @@ -789,3 +789,41 @@ PML::Restart (const std::string& dir) VisMF::Read(*pml_B_cp[2], dir+"_Bz_cp"); } } + +#ifdef WARPX_USE_PSATD +void +PML::PushPSATD() { + SpectralSolver& solver = *(spectral_solver_fp); + + using Idx = SpectralPMLIndex; + + // Perform forward Fourier transform + solver.ForwardTransform(*pml_E_fp[0], Idx::Exy, 0); + solver.ForwardTransform(*pml_E_fp[0], Idx::Exz, 1); + solver.ForwardTransform(*pml_E_fp[1], Idx::Eyx, 0); + solver.ForwardTransform(*pml_E_fp[1], Idx::Eyz, 1); + solver.ForwardTransform(*pml_E_fp[2], Idx::Ezx, 0); + solver.ForwardTransform(*pml_E_fp[2], Idx::Ezy, 1); + solver.ForwardTransform(*pml_B_fp[0], Idx::Bxy, 0); + solver.ForwardTransform(*pml_B_fp[0], Idx::Bxz, 1); + solver.ForwardTransform(*pml_B_fp[1], Idx::Byx, 0); + solver.ForwardTransform(*pml_B_fp[1], Idx::Byz, 1); + solver.ForwardTransform(*pml_B_fp[2], Idx::Bzx, 0); + solver.ForwardTransform(*pml_B_fp[2], Idx::Bzy, 1); + // Advance fields in spectral space + solver.pushSpectralFields(); + // Perform backward Fourier Transform + solver.BackwardTransform(*pml_E_fp[0], Idx::Exy, 0); + solver.BackwardTransform(*pml_E_fp[0], Idx::Exz, 1); + solver.BackwardTransform(*pml_E_fp[1], Idx::Eyx, 0); + solver.BackwardTransform(*pml_E_fp[1], Idx::Eyz, 1); + solver.BackwardTransform(*pml_E_fp[2], Idx::Ezx, 0); + solver.BackwardTransform(*pml_E_fp[2], Idx::Ezy, 1); + solver.BackwardTransform(*pml_B_fp[0], Idx::Bxy, 0); + solver.BackwardTransform(*pml_B_fp[0], Idx::Bxz, 1); + solver.BackwardTransform(*pml_B_fp[1], Idx::Byx, 0); + solver.BackwardTransform(*pml_B_fp[1], Idx::Byz, 1); + solver.BackwardTransform(*pml_B_fp[2], Idx::Bzx, 0); + solver.BackwardTransform(*pml_B_fp[2], Idx::Bzy, 1); +} +#endif diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index 2047a569d..5f0e1d7dd 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -367,17 +367,24 @@ WarpX::PushPSATD (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(dt[lev] == a_dt, "dt must be consistent"); + + // Evolve the fields in the regular boxes if (fft_hybrid_mpi_decomposition){ PushPSATD_hybridFFT(lev, a_dt); } else { PushPSATD_localFFT(lev, a_dt); } + + // Evolve the fields in the PML boxes + if (do_pml) { + pml[lev]->PushPSATD(); + } } } void WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) { - auto& solver = *spectral_solver_fp[lev]; + SpectralSolver& solver = *spectral_solver_fp[lev]; // Perform forward Fourier transform solver.ForwardTransform(*Efield_fp[lev][0], SpectralFieldIndex::Ex); -- cgit v1.2.3 From f13fdc2ac4d1784b5dd6bf4f1172938720c3683f Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Sun, 5 May 2019 16:43:30 -0700 Subject: Correct boxarray --- Source/BoundaryConditions/PML.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 77693915e..30db4d363 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -395,7 +395,10 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, #ifdef WARPX_USE_PSATD const bool in_pml = true; // Tells spectral solver to use split-PML equations const RealVect dx{AMREX_D_DECL(geom->CellSize(0), geom->CellSize(1), geom->CellSize(2))}; - spectral_solver_fp.reset( new SpectralSolver( ba, dm, + // Get the cell-centered box, with guard cells + BoxArray realspace_ba = ba; // Copy box + realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells + spectral_solver_fp.reset( new SpectralSolver( realspace_ba, dm, nox_fft, noy_fft, noz_fft, do_nodal, dx, dt, in_pml ) ); #endif @@ -437,7 +440,10 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, #ifdef WARPX_USE_PSATD const bool in_pml = true; // Tells spectral solver to use split-PML equations const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))}; - spectral_solver_fp.reset( new SpectralSolver( cba, cdm, + // Get the cell-centered box, with guard cells + BoxArray realspace_ba = cba; // Copy box + realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells + spectral_solver_cp.reset( new SpectralSolver( realspace_ba, cdm, nox_fft, noy_fft, noz_fft, do_nodal, cdx, dt, in_pml ) ); #endif } -- cgit v1.2.3 From e53aaff95bdbf1bf5573f775e84640c212b3d121 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Thu, 25 Jul 2019 16:24:12 -0700 Subject: Fix bug in spectral PML --- Source/BoundaryConditions/PML.cpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 5bcec5065..6c00c1d93 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -804,16 +804,20 @@ PML::PushPSATD() { using Idx = SpectralPMLIndex; // Perform forward Fourier transform + // Note: the correspondance between the spectral PML index + // (Exy, Ezx, etc.) and the component (0 or 1) of the + // MultiFabs (e.g. pml_E_fp) is dictated by the + // function that damps the PML solver.ForwardTransform(*pml_E_fp[0], Idx::Exy, 0); solver.ForwardTransform(*pml_E_fp[0], Idx::Exz, 1); - solver.ForwardTransform(*pml_E_fp[1], Idx::Eyx, 0); - solver.ForwardTransform(*pml_E_fp[1], Idx::Eyz, 1); + solver.ForwardTransform(*pml_E_fp[1], Idx::Eyz, 0); + solver.ForwardTransform(*pml_E_fp[1], Idx::Eyx, 1); solver.ForwardTransform(*pml_E_fp[2], Idx::Ezx, 0); solver.ForwardTransform(*pml_E_fp[2], Idx::Ezy, 1); solver.ForwardTransform(*pml_B_fp[0], Idx::Bxy, 0); solver.ForwardTransform(*pml_B_fp[0], Idx::Bxz, 1); - solver.ForwardTransform(*pml_B_fp[1], Idx::Byx, 0); - solver.ForwardTransform(*pml_B_fp[1], Idx::Byz, 1); + solver.ForwardTransform(*pml_B_fp[1], Idx::Byz, 0); + solver.ForwardTransform(*pml_B_fp[1], Idx::Byx, 1); solver.ForwardTransform(*pml_B_fp[2], Idx::Bzx, 0); solver.ForwardTransform(*pml_B_fp[2], Idx::Bzy, 1); // Advance fields in spectral space @@ -821,14 +825,14 @@ PML::PushPSATD() { // Perform backward Fourier Transform solver.BackwardTransform(*pml_E_fp[0], Idx::Exy, 0); solver.BackwardTransform(*pml_E_fp[0], Idx::Exz, 1); - solver.BackwardTransform(*pml_E_fp[1], Idx::Eyx, 0); - solver.BackwardTransform(*pml_E_fp[1], Idx::Eyz, 1); + solver.BackwardTransform(*pml_E_fp[1], Idx::Eyz, 0); + solver.BackwardTransform(*pml_E_fp[1], Idx::Eyx, 1); solver.BackwardTransform(*pml_E_fp[2], Idx::Ezx, 0); solver.BackwardTransform(*pml_E_fp[2], Idx::Ezy, 1); solver.BackwardTransform(*pml_B_fp[0], Idx::Bxy, 0); solver.BackwardTransform(*pml_B_fp[0], Idx::Bxz, 1); - solver.BackwardTransform(*pml_B_fp[1], Idx::Byx, 0); - solver.BackwardTransform(*pml_B_fp[1], Idx::Byz, 1); + solver.BackwardTransform(*pml_B_fp[1], Idx::Byz, 0); + solver.BackwardTransform(*pml_B_fp[1], Idx::Byx, 1); solver.BackwardTransform(*pml_B_fp[2], Idx::Bzx, 0); solver.BackwardTransform(*pml_B_fp[2], Idx::Bzy, 1); } -- cgit v1.2.3 From b7b507d6f521e3bd31bf6ac4f86af4ea12471231 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 30 Jul 2019 10:12:57 -0700 Subject: Put PML spectral-specific arguments in precompiler directives --- Source/BoundaryConditions/PML.H | 6 ++++-- Source/BoundaryConditions/PML.cpp | 6 ++++-- Source/Initialization/WarpXInitData.cpp | 16 +++++++++------- Source/WarpX.H | 13 +++++++------ 4 files changed, 24 insertions(+), 17 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 469805843..65cf73bbf 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -97,8 +97,10 @@ class PML public: PML (const amrex::BoxArray& ba, const amrex::DistributionMapping& dm, const amrex::Geometry* geom, const amrex::Geometry* cgeom, - int ncell, int delta, int ref_ratio, amrex::Real dt, - int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int ncell, int delta, int ref_ratio, +#ifdef WARPX_USE_PSATD + amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, +#endif int do_dive_cleaning, int do_moving_window); void ComputePMLFactors (amrex::Real dt); diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 6c00c1d93..b90d720e8 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -329,8 +329,10 @@ MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, const Geometry* geom, const Geometry* cgeom, - int ncell, int delta, int ref_ratio, Real dt, - int nox_fft, int noy_fft, int noz_fft, bool do_nodal, + int ncell, int delta, int ref_ratio, +#ifdef WARPX_USE_PSATD + Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, +#endif int do_dive_cleaning, int do_moving_window) : m_geom(geom), m_cgeom(cgeom) diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 6f27787e8..332acb473 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -1,6 +1,4 @@ -#include - #include #include @@ -134,15 +132,19 @@ WarpX::InitPML () if (do_pml) { pml[0].reset(new PML(boxArray(0), DistributionMap(0), &Geom(0), nullptr, - pml_ncell, pml_delta, 0, dt[0], - nox_fft, noy_fft, noz_fft, do_nodal, + pml_ncell, pml_delta, 0, +#ifdef WARPX_USE_PSATD + dt[0], nox_fft, noy_fft, noz_fft, do_nodal, +#endif do_dive_cleaning, do_moving_window)); for (int lev = 1; lev <= finest_level; ++lev) { pml[lev].reset(new PML(boxArray(lev), DistributionMap(lev), &Geom(lev), &Geom(lev-1), - pml_ncell, pml_delta, refRatio(lev-1)[0], dt[lev], - nox_fft, noy_fft, noz_fft, do_nodal, + pml_ncell, pml_delta, refRatio(lev-1)[0], +#ifdef WARPX_USE_PSATD + dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, +#endif do_dive_cleaning, do_moving_window)); } } @@ -325,7 +327,7 @@ WarpX::InitLevelData (int lev, Real time) void WarpX::InitLevelDataFFT (int lev, Real time) { - + Efield_fp_fft[lev][0]->setVal(0.0); Efield_fp_fft[lev][1]->setVal(0.0); Efield_fp_fft[lev][2]->setVal(0.0); diff --git a/Source/WarpX.H b/Source/WarpX.H index e9fe642a9..bf10b2a69 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -152,12 +152,12 @@ public: BilinearFilter bilinear_filter; amrex::Vector< std::unique_ptr > nci_godfrey_filter_exeybz; amrex::Vector< std::unique_ptr > nci_godfrey_filter_bxbyez; - + static int num_mirrors; amrex::Vector mirror_z; amrex::Vector mirror_z_width; amrex::Vector mirror_z_npoints; - + void applyMirrors(amrex::Real time); void ComputeDt (); @@ -499,7 +499,7 @@ private: int warpx_do_continuous_injection = 0; int num_injected_species = -1; amrex::Vector injected_plasma_species; - + int do_electrostatic = 0; int n_buffer = 4; amrex::Real const_dt = 0.5e-11; @@ -594,9 +594,6 @@ private: amrex::Vector< std::unique_ptr > rho_cp_fft; #endif -int nox_fft = 16; -int noy_fft = 16; -int noz_fft = 16; #ifdef WARPX_USE_PSATD private: void EvolvePSATD (int numsteps); @@ -606,6 +603,10 @@ private: 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> spectral_solver_fp; amrex::Vector> spectral_solver_cp; #endif -- cgit v1.2.3 From 551bb104a69585cf7560a4aa5901c28f25d231f7 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 30 Jul 2019 10:58:28 -0700 Subject: Apply spectral PML solver to both coarse and fine patch --- Source/BoundaryConditions/PML.H | 6 ++++ Source/BoundaryConditions/PML.cpp | 66 +++++++++++++++++++++++---------------- 2 files changed, 45 insertions(+), 27 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 65cf73bbf..b04938cd8 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -177,4 +177,10 @@ private: static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; +#ifdef WARPX_USE_PSATD +void PushPMLPSATDSinglePatch( SpectralSolver& solver, + std::array,3>& pml_E, + std::array,3>& pml_B ); +#endif + #endif diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index b90d720e8..45a7e360a 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -800,42 +800,54 @@ PML::Restart (const std::string& dir) #ifdef WARPX_USE_PSATD void -PML::PushPSATD() { - SpectralSolver& solver = *(spectral_solver_fp); +PML::PushPSATD () { + + // Update the fields on the fine and coarse patch + PushPMLPSATDSinglePatch( *spectral_solver_fp, pml_E_fp, pml_B_fp ); + if (spectral_solver_cp) { + PushPMLPSATDSinglePatch( *spectral_solver_cp, pml_E_cp, pml_B_cp ); + } +} + +void +PushPMLPSATDSinglePatch ( + SpectralSolver& solver, + std::array,3>& pml_E, + std::array,3>& pml_B ) { using Idx = SpectralPMLIndex; // Perform forward Fourier transform // Note: the correspondance between the spectral PML index // (Exy, Ezx, etc.) and the component (0 or 1) of the - // MultiFabs (e.g. pml_E_fp) is dictated by the + // MultiFabs (e.g. pml_E) is dictated by the // function that damps the PML - solver.ForwardTransform(*pml_E_fp[0], Idx::Exy, 0); - solver.ForwardTransform(*pml_E_fp[0], Idx::Exz, 1); - solver.ForwardTransform(*pml_E_fp[1], Idx::Eyz, 0); - solver.ForwardTransform(*pml_E_fp[1], Idx::Eyx, 1); - solver.ForwardTransform(*pml_E_fp[2], Idx::Ezx, 0); - solver.ForwardTransform(*pml_E_fp[2], Idx::Ezy, 1); - solver.ForwardTransform(*pml_B_fp[0], Idx::Bxy, 0); - solver.ForwardTransform(*pml_B_fp[0], Idx::Bxz, 1); - solver.ForwardTransform(*pml_B_fp[1], Idx::Byz, 0); - solver.ForwardTransform(*pml_B_fp[1], Idx::Byx, 1); - solver.ForwardTransform(*pml_B_fp[2], Idx::Bzx, 0); - solver.ForwardTransform(*pml_B_fp[2], Idx::Bzy, 1); + solver.ForwardTransform(*pml_E[0], Idx::Exy, 0); + solver.ForwardTransform(*pml_E[0], Idx::Exz, 1); + solver.ForwardTransform(*pml_E[1], Idx::Eyz, 0); + solver.ForwardTransform(*pml_E[1], Idx::Eyx, 1); + solver.ForwardTransform(*pml_E[2], Idx::Ezx, 0); + solver.ForwardTransform(*pml_E[2], Idx::Ezy, 1); + solver.ForwardTransform(*pml_B[0], Idx::Bxy, 0); + solver.ForwardTransform(*pml_B[0], Idx::Bxz, 1); + solver.ForwardTransform(*pml_B[1], Idx::Byz, 0); + solver.ForwardTransform(*pml_B[1], Idx::Byx, 1); + solver.ForwardTransform(*pml_B[2], Idx::Bzx, 0); + solver.ForwardTransform(*pml_B[2], Idx::Bzy, 1); // Advance fields in spectral space solver.pushSpectralFields(); // Perform backward Fourier Transform - solver.BackwardTransform(*pml_E_fp[0], Idx::Exy, 0); - solver.BackwardTransform(*pml_E_fp[0], Idx::Exz, 1); - solver.BackwardTransform(*pml_E_fp[1], Idx::Eyz, 0); - solver.BackwardTransform(*pml_E_fp[1], Idx::Eyx, 1); - solver.BackwardTransform(*pml_E_fp[2], Idx::Ezx, 0); - solver.BackwardTransform(*pml_E_fp[2], Idx::Ezy, 1); - solver.BackwardTransform(*pml_B_fp[0], Idx::Bxy, 0); - solver.BackwardTransform(*pml_B_fp[0], Idx::Bxz, 1); - solver.BackwardTransform(*pml_B_fp[1], Idx::Byz, 0); - solver.BackwardTransform(*pml_B_fp[1], Idx::Byx, 1); - solver.BackwardTransform(*pml_B_fp[2], Idx::Bzx, 0); - solver.BackwardTransform(*pml_B_fp[2], Idx::Bzy, 1); + solver.BackwardTransform(*pml_E[0], Idx::Exy, 0); + solver.BackwardTransform(*pml_E[0], Idx::Exz, 1); + solver.BackwardTransform(*pml_E[1], Idx::Eyz, 0); + solver.BackwardTransform(*pml_E[1], Idx::Eyx, 1); + solver.BackwardTransform(*pml_E[2], Idx::Ezx, 0); + solver.BackwardTransform(*pml_E[2], Idx::Ezy, 1); + solver.BackwardTransform(*pml_B[0], Idx::Bxy, 0); + solver.BackwardTransform(*pml_B[0], Idx::Bxz, 1); + solver.BackwardTransform(*pml_B[1], Idx::Byz, 0); + solver.BackwardTransform(*pml_B[1], Idx::Byx, 1); + solver.BackwardTransform(*pml_B[2], Idx::Bzx, 0); + solver.BackwardTransform(*pml_B[2], Idx::Bzy, 1); } #endif -- cgit v1.2.3 From 90c7a587a87e4411434c0dc6433bc6c5703fe521 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Tue, 6 Aug 2019 10:43:48 -0700 Subject: Update naming in PML PSATD solver --- Source/BoundaryConditions/PML.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 45a7e360a..cd723f889 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -443,9 +443,9 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, const bool in_pml = true; // Tells spectral solver to use split-PML equations const RealVect cdx{AMREX_D_DECL(cgeom->CellSize(0), cgeom->CellSize(1), cgeom->CellSize(2))}; // Get the cell-centered box, with guard cells - BoxArray realspace_ba = cba; // Copy box - realspace_ba.enclosedCells().grow(nge); // cell-centered + guard cells - spectral_solver_cp.reset( new SpectralSolver( realspace_ba, cdm, + BoxArray realspace_cba = cba; // Copy box + realspace_cba.enclosedCells().grow(nge); // cell-centered + guard cells + spectral_solver_cp.reset( new SpectralSolver( realspace_cba, cdm, nox_fft, noy_fft, noz_fft, do_nodal, cdx, dt, in_pml ) ); #endif } -- cgit v1.2.3 From 4150e6efc4ee7c0d6493c6f2464135623b075773 Mon Sep 17 00:00:00 2001 From: Remi Lehe Date: Wed, 7 Aug 2019 15:41:08 -0700 Subject: Clean-up if conditions in PML --- Source/BoundaryConditions/PML.cpp | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index cd723f889..96bc08af9 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -258,14 +258,7 @@ SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) { for (int i = 0, N = sigma_star[idim].size(); i < N; ++i) { - if (sigma_star[idim][i] == 0.0) - { - sigma_star_fac[idim][i] = 1.0; - } - else - { - sigma_star_fac[idim][i] = std::exp(-sigma_star[idim][i]*dt); - } + sigma_star_fac[idim][i] = std::exp(-sigma_star[idim][i]*dt); } } } @@ -277,14 +270,7 @@ SigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) { for (int i = 0, N = sigma[idim].size(); i < N; ++i) { - if (sigma[idim][i] == 0.0) - { - sigma_fac[idim][i] = 1.0; - } - else - { - sigma_fac[idim][i] = std::exp(-sigma[idim][i]*dt); - } + sigma_fac[idim][i] = std::exp(-sigma[idim][i]*dt); } } } -- cgit v1.2.3 From 23f57db568dea4dc6545482e4846b35102558fb4 Mon Sep 17 00:00:00 2001 From: ablelly Date: Thu, 8 Aug 2019 01:26:36 +0200 Subject: Made all changes in function calls to include do_pml_Lo, do_pml_Hi --- Source/BoundaryConditions/PML.H | 8 ++++++-- Source/BoundaryConditions/PML.cpp | 6 ++++-- Source/Initialization/WarpXInitData.cpp | 3 ++- 3 files changed, 12 insertions(+), 5 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index b04938cd8..7ded63fc3 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -101,7 +101,9 @@ public: #ifdef WARPX_USE_PSATD amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, #endif - int do_dive_cleaning, int do_moving_window); + int do_dive_cleaning, int do_moving_window, + const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), + const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); void ComputePMLFactors (amrex::Real dt); @@ -172,7 +174,9 @@ private: #endif static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, - const amrex::BoxArray& grid_ba, int ncell); + const amrex::BoxArray& grid_ba, int ncell, + const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), + const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 96bc08af9..e444e3b5b 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -319,7 +319,8 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, #ifdef WARPX_USE_PSATD Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, #endif - int do_dive_cleaning, int do_moving_window) + int do_dive_cleaning, int do_moving_window, + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) : m_geom(geom), m_cgeom(cgeom) { @@ -438,7 +439,8 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, } BoxArray -PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell) +PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell, + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi) { Box domain = geom.Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 332acb473..f93c3705e 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -145,7 +145,8 @@ WarpX::InitPML () #ifdef WARPX_USE_PSATD dt[lev], nox_fft, noy_fft, noz_fft, do_nodal, #endif - do_dive_cleaning, do_moving_window)); + do_dive_cleaning, do_moving_window, + do_pml_Lo, do_pml_Hi)); } } } -- cgit v1.2.3 From 1d8f4491b66d4d58f45986da119ba23f6ad1a85e Mon Sep 17 00:00:00 2001 From: ablelly Date: Thu, 8 Aug 2019 01:36:43 +0200 Subject: Added the specificity do_pml_Lo, do_pml_Hi --- Source/BoundaryConditions/PML.H | 7 +++---- Source/BoundaryConditions/PML.cpp | 11 ++++++++--- Source/Initialization/WarpXInitData.cpp | 3 ++- 3 files changed, 13 insertions(+), 8 deletions(-) (limited to 'Source/BoundaryConditions/PML.cpp') diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 7ded63fc3..b34cbe88b 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -102,8 +102,7 @@ public: amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, bool do_nodal, #endif int do_dive_cleaning, int do_moving_window, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const amrex::IntVect do_pml_Lo, const amrex::IntVect do_pml_Hi); void ComputePMLFactors (amrex::Real dt); @@ -175,8 +174,8 @@ private: static amrex::BoxArray MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell, - const amrex::IntVect do_pml_Lo = amrex::IntVect::TheUnitVector(), - const amrex::IntVect do_pml_Hi = amrex::IntVect::TheUnitVector()); + const amrex::IntVect do_pml_Lo, + const amrex::IntVect do_pml_Hi); static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index e444e3b5b..21d348482 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -324,7 +324,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, : m_geom(geom), m_cgeom(cgeom) { - const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell); + const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell, do_pml_Lo, do_pml_Hi); if (ba.size() == 0) { m_ok = false; return; @@ -400,7 +400,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, BoxArray grid_cba = grid_ba; grid_cba.coarsen(ref_ratio); - const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell); + const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell, do_pml_Lo, do_pml_Hi); DistributionMapping cdm{cba}; @@ -445,7 +445,12 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, Box domain = geom.Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { if ( ! geom.isPeriodic(idim) ) { - domain.grow(idim, ncell); + if (do_pml_Lo[idim]){ + domain.growLo(idim, ncell); + } + if (do_pml_Hi[idim]){ + domain.growHi(idim, ncell); + } } } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index f93c3705e..7ca8e14c4 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -136,7 +136,8 @@ WarpX::InitPML () #ifdef WARPX_USE_PSATD dt[0], nox_fft, noy_fft, noz_fft, do_nodal, #endif - do_dive_cleaning, do_moving_window)); + do_dive_cleaning, do_moving_window, + do_pml_Lo, do_pml_Hi)); for (int lev = 1; lev <= finest_level; ++lev) { pml[lev].reset(new PML(boxArray(lev), DistributionMap(lev), -- cgit v1.2.3