diff options
| -rw-r--r-- | .gitignore | 2 | ||||
| -rw-r--r-- | Examples/Tests/PML/inputs2d | 2 | ||||
| -rw-r--r-- | Source/BoundaryConditions/PML.H | 52 | ||||
| -rw-r--r-- | Source/BoundaryConditions/PML.cpp | 411 | ||||
| -rw-r--r-- | Source/BoundaryConditions/PML_Cplx.cpp | 958 | ||||
| -rw-r--r-- | Source/BoundaryConditions/PML_routines.F90 | 482 | ||||
| -rw-r--r-- | Source/BoundaryConditions/WarpXEvolvePML.cpp | 53 | ||||
| -rw-r--r-- | Source/Evolve/WarpXEvolveEM.cpp | 80 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 32 | ||||
| -rw-r--r-- | Source/FortranInterface/WarpX_f.H | 42 | ||||
| -rw-r--r-- | Source/Initialization/WarpXInitData.cpp | 16 | ||||
| -rw-r--r-- | Source/WarpX.H | 12 | ||||
| -rw-r--r-- | Source/WarpX.cpp | 2 |
13 files changed, 2015 insertions, 129 deletions
diff --git a/.gitignore b/.gitignore index 02bd9eafa..0a57eb54c 100644 --- a/.gitignore +++ b/.gitignore @@ -9,4 +9,4 @@ AMReX_buildInfo.cpp Backtrace.* *.vth *.sw[opqrs] - +Bin/ diff --git a/Examples/Tests/PML/inputs2d b/Examples/Tests/PML/inputs2d index 5b936a333..dc1f136ef 100644 --- a/Examples/Tests/PML/inputs2d +++ b/Examples/Tests/PML/inputs2d @@ -2,7 +2,7 @@ max_step = 300 # number of grid points -amr.n_cell = 128 512 +amr.n_cell = 128 256 # Maximum allowable size of each subdomain in the problem domain; # this is used to decompose the domain for parallel calculations. diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index 0cf367284..fd8fe12ec 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -16,6 +16,22 @@ (x).sigma_star_fac[1].data(), (x).sigma_star_fac[1].m_lo, (x).sigma_star_fac[1].m_hi, \ (x).sigma_star_fac[2].data(), (x).sigma_star_fac[2].m_lo, (x).sigma_star_fac[2].m_hi +#define WRPX_PML_SIGMAJ_TO_FORTRAN(x) \ + (x).sigma_cum_fac[0].data(), (x).sigma_cum_fac[0].m_lo, (x).sigma_cum_fac[0].m_hi, \ + (x).sigma_cum_fac[1].data(), (x).sigma_cum_fac[1].m_lo, (x).sigma_cum_fac[1].m_hi, \ + (x).sigma_cum_fac[2].data(), (x).sigma_cum_fac[2].m_lo, (x).sigma_cum_fac[2].m_hi, \ + (x).sigma_star_cum_fac[0].data(), (x).sigma_star_cum_fac[0].m_lo, (x).sigma_star_cum_fac[0].m_hi, \ + (x).sigma_star_cum_fac[1].data(), (x).sigma_star_cum_fac[1].m_lo, (x).sigma_star_cum_fac[1].m_hi, \ + (x).sigma_star_cum_fac[2].data(), (x).sigma_star_cum_fac[2].m_lo, (x).sigma_star_cum_fac[2].m_hi + +#define WRPX_PML_SIGMACOEFF_TO_FORTRAN(x) \ + (x).sigma[0].dataPtr(), &((x).sigma[0].m_lo), &((x).sigma[0].m_hi), \ + (x).sigma[1].dataPtr(), &((x).sigma[1].m_lo), &((x).sigma[1].m_hi), \ + (x).sigma[2].dataPtr(), &((x).sigma[2].m_lo), &((x).sigma[2].m_hi), \ + (x).sigma_star[0].dataPtr(), &((x).sigma_star[0].m_lo), &((x).sigma_star[0].m_hi), \ + (x).sigma_star[1].dataPtr(), &((x).sigma_star[1].m_lo), &((x).sigma_star[1].m_hi), \ + (x).sigma_star[2].dataPtr(), &((x).sigma_star[2].m_lo), &((x).sigma_star[2].m_hi) + #else #define WRPX_PML_TO_FORTRAN(x) \ @@ -24,6 +40,18 @@ (x).sigma_star_fac[0].data(), (x).sigma_star_fac[0].m_lo, (x).sigma_star_fac[0].m_hi, \ (x).sigma_star_fac[1].data(), (x).sigma_star_fac[1].m_lo, (x).sigma_star_fac[1].m_hi +#define WRPX_PML_SIGMAJ_TO_FORTRAN(x) \ + (x).sigma_cum_fac[0].data(), (x).sigma_cum_fac[0].m_lo, (x).sigma_cum_fac[0].m_hi, \ + (x).sigma_cum_fac[1].data(), (x).sigma_cum_fac[1].m_lo, (x).sigma_cum_fac[1].m_hi, \ + (x).sigma_star_cum_fac[0].data(), (x).sigma_star_cum_fac[0].m_lo, (x).sigma_star_cum_fac[0].m_hi, \ + (x).sigma_star_cum_fac[1].data(), (x).sigma_star_cum_fac[1].m_lo, (x).sigma_star_cum_fac[1].m_hi + +#define WRPX_PML_SIGMACOEFF_TO_FORTRAN(x) \ + (x).sigma[0].data(), (x).sigma[0].m_lo, (x).sigma[0].m_hi, \ + (x).sigma[1].data(), (x).sigma[1].m_lo, (x).sigma[1].m_hi, \ + (x).sigma_star[0].data(), (x).sigma_star[0].m_lo, (x).sigma_star[0].m_hi, \ + (x).sigma_star[1].data(), (x).sigma_star[1].m_lo, (x).sigma_star[1].m_hi + #endif struct Sigma : amrex::Vector<amrex::Real> @@ -44,9 +72,16 @@ struct SigmaBox using SigmaVect = std::array<Sigma,AMREX_SPACEDIM>; SigmaVect sigma; // sigma/epsilon + SigmaVect sigma_cum; // cumsum(sigma)/(c*epsilon) SigmaVect sigma_star; // sigma_star/mu + SigmaVect sigma_star_cum; SigmaVect sigma_fac; + SigmaVect sigma_cum_fac; SigmaVect sigma_star_fac; + SigmaVect sigma_star_cum_fac; + + std::string pml_type; // if pml normal to x : if pml Hi : pml_type={1,0,0}; if pml Lo: pml_type={-1,0,0} + std::array<int,5> pml_type_array; }; namespace amrex { @@ -93,14 +128,16 @@ 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, int do_dive_cleaning, int do_moving_window, int pml_has_particles); void ComputePMLFactors (amrex::Real dt); std::array<amrex::MultiFab*,3> GetE_fp (); std::array<amrex::MultiFab*,3> GetB_fp (); + std::array<amrex::MultiFab*,3> Getj_fp (); std::array<amrex::MultiFab*,3> GetE_cp (); std::array<amrex::MultiFab*,3> GetB_cp (); + std::array<amrex::MultiFab*,3> Getj_cp (); amrex::MultiFab* GetF_fp (); amrex::MultiFab* GetF_cp (); @@ -115,10 +152,19 @@ 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 CopyJinPMLs (const std::array<amrex::MultiFab*,3>& j_fp, + const std::array<amrex::MultiFab*,3>& j_cp); + void CopyJinReg (const std::array<amrex::MultiFab*,3>& j_fp, + const std::array<amrex::MultiFab*,3>& j_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 CopyJinPMLs (PatchType patch_type, + const std::array<amrex::MultiFab*,3>& jp); + void CopyJinReg (PatchType patch_type, + const std::array<amrex::MultiFab*,3>& jp); void ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp); void ExchangeF (PatchType patch_type, amrex::MultiFab* Fp); @@ -144,9 +190,11 @@ private: std::array<std::unique_ptr<amrex::MultiFab>,3> pml_E_fp; std::array<std::unique_ptr<amrex::MultiFab>,3> pml_B_fp; + std::array<std::unique_ptr<amrex::MultiFab>,3> pml_j_fp; std::array<std::unique_ptr<amrex::MultiFab>,3> pml_E_cp; std::array<std::unique_ptr<amrex::MultiFab>,3> pml_B_cp; + std::array<std::unique_ptr<amrex::MultiFab>,3> pml_j_cp; std::unique_ptr<amrex::MultiFab> pml_F_fp; std::unique_ptr<amrex::MultiFab> pml_F_cp; @@ -158,6 +206,8 @@ private: const amrex::BoxArray& grid_ba, int ncell); static void Exchange (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); + static void CopyRegInPMLs (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); + static void CopyPMLsInReg (amrex::MultiFab& pml, amrex::MultiFab& reg, const amrex::Geometry& geom); }; #endif diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index f780f335c..b01a95f10 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -1,4 +1,3 @@ - #include <PML.H> #include <WarpX.H> #include <WarpXConst.H> @@ -16,59 +15,78 @@ using namespace amrex; namespace { - static void FillLo (int idim, Sigma& sigma, Sigma& sigma_star, + static void FillLo (int idim, Sigma& sigma, Sigma& sigma_cum, Sigma& sigma_star, Sigma& sigma_star_cum, const Box& overlap, const Box& grid, Real fac) { int glo = grid.smallEnd(idim); int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; + int shi = sigma.m_hi; int sslo = sigma_star.m_lo; + Real x = 10.0; + Real theta = 10.0; + Real coeff_damp = 1.; //std::sin(theta*MathConst::pi/180.); + for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast<Real>(glo-i); sigma[i-slo] = fac*(offset*offset); + sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); } + for (int i = olo; i <= ohi; ++i) { Real offset = static_cast<Real>(glo-i) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); + sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); } } - static void FillHi (int idim, Sigma& sigma, Sigma& sigma_star, + static void FillHi (int idim, Sigma& sigma, Sigma& sigma_cum, Sigma& sigma_star, Sigma& sigma_star_cum, const Box& overlap, const Box& grid, Real fac) { int ghi = grid.bigEnd(idim); int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; + int shi = sigma.m_hi; int sslo = sigma_star.m_lo; + Real x = 10.0; + Real theta = 10.0; + Real coeff_damp = 1.;//std::sin(theta*MathConst::pi/180.); for (int i = olo; i <= ohi+1; ++i) { Real offset = static_cast<Real>(i-ghi-1); sigma[i-slo] = fac*(offset*offset); + sigma_cum[i-slo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); } for (int i = olo; i <= ohi; ++i) { Real offset = static_cast<Real>(i-ghi) - 0.5; sigma_star[i-sslo] = fac*(offset*offset); + sigma_star_cum[i-sslo] = coeff_damp*(fac*(offset*offset*offset)/3.)/(PhysConst::c*x/std::sqrt(1+x*x)); } } - static void FillZero (int idim, Sigma& sigma, Sigma& sigma_star, const Box& overlap) + static void FillZero (int idim, Sigma& sigma, Sigma& sigma_cum, Sigma& sigma_star, Sigma& sigma_star_cum, const Box& overlap) { int olo = overlap.smallEnd(idim); int ohi = overlap.bigEnd(idim); int slo = sigma.m_lo; int sslo = sigma_star.m_lo; std::fill(sigma.begin()+(olo-slo), sigma.begin()+(ohi+2-slo), 0.0); + std::fill(sigma_cum.begin()+(olo-slo), sigma_cum.begin()+(ohi+2-slo), 0.0); std::fill(sigma_star.begin()+(olo-sslo), sigma_star.begin()+(ohi+1-sslo), 0.0); + std::fill(sigma_star_cum.begin()+(olo-sslo), sigma_star_cum.begin()+(ohi+1-sslo), 0.0); } } SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta) { + amrex::Print()<< "===== BUILD SigmaBox ====="<<std::endl; + // amrex::Print()<<"box = ["<<box.smallEnd()[0]<<", "<<box.smallEnd()[1]<<", "<<box.bigEnd()[0]<<", "<<box.bigEnd()[1]<<"]"<<std::endl; + amrex::Print()<<"box = ["<<box.smallEnd()[0]<<", "<<box.smallEnd()[1]<<", "<<box.smallEnd()[2]<<", "<<box.bigEnd()[0]<<", "<<box.bigEnd()[1]<<", "<<box.bigEnd()[2]<<"]"<<std::endl; BL_ASSERT(box.cellCentered()); const IntVect& sz = box.size(); @@ -78,20 +96,35 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { sigma [idim].resize(sz[idim]+1); - sigma_star [idim].resize(sz[idim] ); + sigma_cum [idim].resize(sz[idim]+1); + sigma_star [idim].resize(sz[idim]); + sigma_star_cum[idim].resize(sz[idim]); sigma_fac [idim].resize(sz[idim]+1); - sigma_star_fac[idim].resize(sz[idim] ); + sigma_cum_fac [idim].resize(sz[idim]+1); + sigma_star_fac[idim].resize(sz[idim]); + sigma_star_cum_fac[idim].resize(sz[idim]); sigma [idim].m_lo = lo[idim]; sigma [idim].m_hi = hi[idim]+1; + sigma_cum [idim].m_lo = lo[idim]; + sigma_cum [idim].m_hi = hi[idim]+1; sigma_star [idim].m_lo = lo[idim]; sigma_star [idim].m_hi = hi[idim]; + sigma_star_cum[idim].m_lo = lo[idim]; + sigma_star_cum[idim].m_hi = hi[idim]; sigma_fac [idim].m_lo = lo[idim]; sigma_fac [idim].m_hi = hi[idim]+1; + sigma_cum_fac [idim].m_lo = lo[idim]; + sigma_cum_fac [idim].m_hi = hi[idim]+1; sigma_star_fac[idim].m_lo = lo[idim]; sigma_star_fac[idim].m_hi = hi[idim]; + sigma_star_cum_fac[idim].m_lo = lo[idim]; + sigma_star_cum_fac[idim].m_hi = hi[idim]; } + pml_type = ""; + pml_type_array = {3, 3, 3, 3, 3}; + Array<Real,AMREX_SPACEDIM> fac; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { fac[idim] = 4.0*PhysConst::c/(dx[idim]*static_cast<Real>(delta*delta)); @@ -114,35 +147,42 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n if (amrex::grow(grid_box, idim, ncell).intersects(box)) { direct_faces.push_back(kv.first); + amrex::Print()<<"direct_faces"<<std::endl; } else if (amrex::grow(grid_box, jdim, ncell).intersects(box)) { side_faces.push_back(kv.first); + amrex::Print()<<"side_faces"<<std::endl; } #if (AMREX_SPACEDIM == 3) else if (amrex::grow(grid_box, kdim, ncell).intersects(box)) { side_faces.push_back(kv.first); + amrex::Print()<<"side_faces 2"<<std::endl; } else if (amrex::grow(amrex::grow(grid_box,idim,ncell), jdim,ncell).intersects(box)) { direct_side_edges.push_back(kv.first); + amrex::Print()<<"direct_side_edges"<<std::endl; } else if (amrex::grow(amrex::grow(grid_box,idim,ncell), kdim,ncell).intersects(box)) { direct_side_edges.push_back(kv.first); + amrex::Print()<<"direct_side_edges 2"<<std::endl; } else if (amrex::grow(amrex::grow(grid_box,jdim,ncell), kdim,ncell).intersects(box)) { side_side_edges.push_back(kv.first); + amrex::Print()<<"side_side_edges"<<std::endl; } #endif else { corners.push_back(kv.first); + amrex::Print()<<"corners"<<std::endl; } } @@ -157,7 +197,20 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n #endif Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + FillLo(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], looverlap, grid_box, fac[idim]); + if (idim == 0){ + pml_type = "crn-++"; // for 2d only use the two first components + pml_type_array[0] = 1; + pml_type_array[1] = 0; + pml_type_array[2] = 1; + pml_type_array[3] = 1; + } + else { + if (pml_type[0] == 'c'){ + pml_type[3+idim] = '-'; + pml_type_array[1+idim] = 0; + } + } } Box hibox = amrex::adjCellHi(grid_box, idim, ncell); @@ -167,7 +220,20 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n #endif Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + FillHi(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], hioverlap, grid_box, fac[idim]); + if (idim == 0){ + pml_type = "crn+++"; // for 2d only use the two first components + pml_type_array[0] = 1; + pml_type_array[1] = 1; + pml_type_array[2] = 1; + pml_type_array[3] = 1; + } + else { + if (pml_type[0] == 'c'){ + pml_type[3+idim] = '+'; + pml_type_array[1+idim] = 1; + } + } } if (!looverlap.ok() && !hioverlap.ok()) { @@ -181,8 +247,13 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& grid_box = grids[gid]; const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_star[idim], overlap); - } else { + FillZero(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], overlap); + if (idim==0){ + pml_type="sed33--"; + pml_type_array[0] = 2; + } + } + else { amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n"); } } @@ -194,13 +265,57 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n Box lobox = amrex::adjCellLo(grid_box, idim, ncell); Box looverlap = lobox.grow(jdim,ncell).grow(kdim,ncell) & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + FillLo(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], looverlap, grid_box, fac[idim]); + if (idim == 0){ + pml_type = "sed03-+"; // for 2d only use the two first components + pml_type_array[0] = 2; + pml_type_array[1] = idim; + pml_type_array[3] = 0; + } + else { + if (pml_type[0] == 's'){ + if (pml_type[3]=='3'){ + pml_type[3]=idim+'0'; + pml_type[5]='-'; + pml_type_array[1] = idim; + pml_type_array[3] = 0; + } + else{ + pml_type[4] = idim+'0'; + pml_type[6] = '-'; + pml_type_array[2] = idim; + pml_type_array[4] = 0; + } + } + } } Box hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox.grow(jdim,ncell).grow(kdim,ncell) & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + FillHi(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], hioverlap, grid_box, fac[idim]); + if (idim == 0){ + pml_type = "sed03++"; // for 2d only use the two first components + pml_type_array[0] = 2; + pml_type_array[1] = idim; + pml_type_array[3] = 1; + } + else { + if (pml_type[0] == 's'){ + if (pml_type[3]=='3'){ + pml_type[3]=idim+'0'; + pml_type[5]='+'; + pml_type_array[1] = idim; + pml_type_array[3] = 1; + } + else{ + pml_type[4] = idim+'0'; + pml_type[6] = '+'; + pml_type_array[2] = idim; + pml_type_array[4] = 1; + } + } + } } if (!looverlap.ok() && !hioverlap.ok()) { @@ -218,7 +333,7 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; #endif if (overlap.ok()) { - FillZero(idim, sigma[idim], sigma_star[idim], overlap); + FillZero(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], overlap); } else { amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n"); } @@ -231,13 +346,23 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell); Box looverlap = lobox & box; if (looverlap.ok()) { - FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + FillLo(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], looverlap, grid_box, fac[idim]); + pml_type = "df0-"; + pml_type_array[0] = 0; + pml_type_array[1] = idim; + pml_type_array[2] = 0; + if (idim > 0){pml_type[2]=idim+'0';} //std::to_string(idim) } const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell); Box hioverlap = hibox & box; if (hioverlap.ok()) { - FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + FillHi(idim, sigma[idim], sigma_cum[idim], sigma_star[idim], sigma_star_cum[idim], hioverlap, grid_box, fac[idim]); + pml_type = "df0+"; + pml_type_array[0] = 0; + pml_type_array[1] = idim; + pml_type_array[2] = 1; + if (idim>0){pml_type[2] = idim+'0';} } if (!looverlap.ok() && !hioverlap.ok()) { @@ -249,8 +374,15 @@ SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int n amrex::Abort("SigmaBox::SigmaBox(): direct_faces.size() > 1, Box gaps not wide enough?\n"); } } + amrex::Print()<<"pml_type = "<<pml_type<<std::endl; + amrex::Print()<<"pml_type_array = ["; + for (int i=0; i<5; i++){ + amrex::Print()<<pml_type_array[i]<<", "; + } + amrex::Print()<<"]"<<std::endl; } + void SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) { @@ -266,6 +398,15 @@ SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) { sigma_star_fac[idim][i] = std::exp(-sigma_star[idim][i]*dt); } + + if (sigma_star_cum[idim][i] == 0.0) + { + sigma_star_cum_fac[idim][i] = 1.0; + } + else + { + sigma_star_cum_fac[idim][i] = std::exp(-sigma_star_cum[idim][i]*dx[idim]); + } } } } @@ -285,6 +426,13 @@ SigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) { sigma_fac[idim][i] = std::exp(-sigma[idim][i]*dt); } + if (sigma_cum[idim][i] == 0.0) + { + sigma_cum_fac[idim][i] = 1.0; + } + else { + sigma_cum_fac[idim][i] = std::exp(-sigma_cum[idim][i]*dx[idim]); + } } } } @@ -329,11 +477,39 @@ 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, int do_dive_cleaning, int do_moving_window, int pml_has_particles) : m_geom(geom), m_cgeom(cgeom) { - const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell); + // BoxList bl_init = BoxList(grid_ba); + // + // amrex::Print() << "========== Printing grid_ba boxes" << std::endl; + // amrex::Print() << "[" << std::endl; + // for (const Box& b: bl_init) { + // amrex::Print() << "[" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "<<b.bigEnd()[0] << ", "<< b.bigEnd()[1] << "]," << std::endl; + // } + // amrex::Print()<< "];" << std::endl; + amrex::Print()<<"===== BUILDING PML ====="<<std::endl; + Box domain0 = geom->Domain(); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( ! geom->isPeriodic(idim) ) { + domain0.grow(idim, -ncell); + } + } + // Box domain0 = amrex::grow(geom->Domain(), -ncell); + amrex::Print() << "[" << domain0.smallEnd()[0]<<", "<< domain0.smallEnd()[1]<<", "<< domain0.smallEnd()[2]<< ", "<<domain0.bigEnd()[0] << ", "<< domain0.bigEnd()[1]<< ", "<< domain0.bigEnd()[2] << "]," << std::endl; + const BoxArray grid_ba_reduced = BoxArray(grid_ba.boxList().intersect(domain0)); + + // BoxList bl_reduced = BoxList(grid_ba_reduced); + // + // amrex::Print() << "========== Printing grid_ba_reduced boxes" << std::endl; + // amrex::Print() << "[" << std::endl; + // for (const Box& b: bl_reduced) { + // amrex::Print() << "[" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "<<b.bigEnd()[0] << ", "<< b.bigEnd()[1] << "]," << std::endl; + // } + // amrex::Print()<< "];" << std::endl; + + const BoxArray& ba = MakeBoxArray(*geom, grid_ba_reduced, ncell); //MakeBoxArray(*geom, grid_ba, ncell); if (ba.size() == 0) { m_ok = false; return; @@ -355,6 +531,7 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_B_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::By_nodal_flag), dm, 2, ngb)); pml_B_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Bz_nodal_flag), dm, 2, ngb)); + pml_E_fp[0]->setVal(0.0); pml_E_fp[1]->setVal(0.0); pml_E_fp[2]->setVal(0.0); @@ -362,13 +539,32 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, pml_B_fp[1]->setVal(0.0); pml_B_fp[2]->setVal(0.0); + pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jx_nodal_flag) + pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jy_nodal_flag) + pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jz_nodal_flag) + pml_j_fp[0]->setVal(0.0); + pml_j_fp[1]->setVal(0.0); + pml_j_fp[2]->setVal(0.0); + + // if (pml_has_particles){ + // pml_j_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::jx_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jx_nodal_flag) + // pml_j_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::jy_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jy_nodal_flag) + // pml_j_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::jz_nodal_flag), dm, 1, ngb)); //convert(ba,WarpX::Jz_nodal_flag) + // pml_j_fp[0]->setVal(0.0); + // pml_j_fp[1]->setVal(0.0); + // pml_j_fp[2]->setVal(0.0); + // amrex::Print() << "PML HAS PARTICLES - fine"<< std::endl; + // + // } + if (do_dive_cleaning) { pml_F_fp.reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()), dm, 3, ngf)); pml_F_fp->setVal(0.0); } - sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); + + sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba_reduced, geom->CellSize(), ncell, delta)); if (cgeom) { @@ -378,7 +574,8 @@ 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 grid_cba_reduced = BoxArray(grid_cba.boxList().intersect(domain0)); + const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba_reduced, ncell); DistributionMapping cdm{cba}; @@ -400,9 +597,27 @@ PML::PML (const BoxArray& grid_ba, const DistributionMapping& grid_dm, { pml_F_cp.reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()), cdm, 3, ngf)); pml_F_cp->setVal(0.0); - } - sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + } + pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb)); + pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb)); + pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb)); + pml_j_cp[0]->setVal(0.0); + pml_j_cp[1]->setVal(0.0); + pml_j_cp[2]->setVal(0.0); + // if (pml_has_particles) + // { + // pml_j_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::jx_nodal_flag), cdm, 1, ngb)); + // pml_j_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::jy_nodal_flag), cdm, 1, ngb)); + // pml_j_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::jz_nodal_flag), cdm, 1, ngb)); + // pml_j_cp[0]->setVal(0.0); + // pml_j_cp[1]->setVal(0.0); + // pml_j_cp[2]->setVal(0.0); + // amrex::Print() << "PML HAS PARTICLES - coarse"<< std::endl; + // } + + // sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba_reduced, cgeom->CellSize(), ncell, delta)); } } @@ -422,8 +637,8 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, { const Box& grid_bx = grid_ba[i]; const IntVect& grid_bx_sz = grid_bx.size(); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, - "Consider using larger amr.blocking_factor"); + // AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, + // "Consider using larger amr.blocking_factor"); Box bx = grid_bx; bx.grow(ncell); @@ -464,6 +679,15 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, BoxArray ba(bl); ba.removeOverlap(false); + // BoxList bl_2 = BoxList(ba); + // + // amrex::Print() << "Printing PML boxes AFTER cleaning" << std::endl; + // amrex::Print() << "[" << std::endl; + // for (const Box& b: bl_2) { + // amrex::Print() << "[" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "<<b.bigEnd()[0] << ", "<< b.bigEnd()[1] << "]," << std::endl; + // } + // amrex::Print()<< "];" << std::endl; + return ba; } @@ -493,6 +717,12 @@ PML::GetB_fp () } std::array<MultiFab*,3> +PML::Getj_fp () +{ + return {pml_j_fp[0].get(), pml_j_fp[1].get(), pml_j_fp[2].get()}; +} + +std::array<MultiFab*,3> PML::GetE_cp () { return {pml_E_cp[0].get(), pml_E_cp[1].get(), pml_E_cp[2].get()}; @@ -504,6 +734,12 @@ PML::GetB_cp () return {pml_B_cp[0].get(), pml_B_cp[1].get(), pml_B_cp[2].get()}; } +std::array<MultiFab*,3> +PML::Getj_cp () +{ + return {pml_j_cp[0].get(), pml_j_cp[1].get(), pml_j_cp[2].get()}; +} + MultiFab* PML::GetF_fp () { @@ -569,6 +805,58 @@ PML::ExchangeE (PatchType patch_type, } void +PML::CopyJinPMLs (PatchType patch_type, + const std::array<amrex::MultiFab*,3>& jp) +{ + if (patch_type == PatchType::fine && pml_j_fp[0] && jp[0]) + { + CopyRegInPMLs(*pml_j_fp[0], *jp[0], *m_geom); + CopyRegInPMLs(*pml_j_fp[1], *jp[1], *m_geom); + CopyRegInPMLs(*pml_j_fp[2], *jp[2], *m_geom); + } + else if (patch_type == PatchType::coarse && pml_j_cp[0] && jp[0]) + { + CopyRegInPMLs(*pml_j_cp[0], *jp[0], *m_cgeom); + CopyRegInPMLs(*pml_j_cp[1], *jp[1], *m_cgeom); + CopyRegInPMLs(*pml_j_cp[2], *jp[2], *m_cgeom); + } +} + +void +PML::CopyJinPMLs (const std::array<amrex::MultiFab*,3>& j_fp, + const std::array<amrex::MultiFab*,3>& j_cp) +{ + CopyJinPMLs(PatchType::fine, j_fp); + CopyJinPMLs(PatchType::coarse, j_cp); +} + +void +PML::CopyJinReg (PatchType patch_type, + const std::array<amrex::MultiFab*,3>& jp) +{ + if (patch_type == PatchType::fine && pml_j_fp[0] && jp[0]) + { + CopyPMLsInReg(*pml_j_fp[0], *jp[0], *m_geom); + CopyPMLsInReg(*pml_j_fp[1], *jp[1], *m_geom); + CopyPMLsInReg(*pml_j_fp[2], *jp[2], *m_geom); + } + else if (patch_type == PatchType::coarse && pml_j_cp[0] && jp[0]) + { + CopyPMLsInReg(*pml_j_cp[0], *jp[0], *m_cgeom); + CopyPMLsInReg(*pml_j_cp[1], *jp[1], *m_cgeom); + CopyPMLsInReg(*pml_j_cp[2], *jp[2], *m_cgeom); + } +} + +void +PML::CopyJinReg (const std::array<amrex::MultiFab*,3>& j_fp, + const std::array<amrex::MultiFab*,3>& j_cp) +{ + CopyJinReg(PatchType::fine, j_fp); + CopyJinReg(PatchType::coarse, j_cp); +} + +void PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp) { ExchangeF(PatchType::fine, F_fp); @@ -588,44 +876,81 @@ PML::ExchangeF (PatchType patch_type, MultiFab* Fp) void PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) { + const IntVect& ngr = reg.nGrowVect(); const IntVect& ngp = pml.nGrowVect(); const int ncp = pml.nComp(); const auto& period = geom.periodicity(); MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), ncp, ngr); + tmpregmf.setVal(0.0, 0, ncp, ngr); + MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), ncp, ngp); + totpmlmf.setVal(0.0, 0, ncp, ngp); + // realise sum of splitted fields inside pml + MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); + if (ncp == 3) { + MultiFab::Add(totpmlmf,pml,2,0,1,0); + } + totpmlmf.setVal(0.0, 1, ncp-1, 0); + reg.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period); if (ngp.max() > 0) // Copy from pml to the ghost cells of regular data { - MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); - MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); - if (ncp == 3) { - MultiFab::Add(totpmlmf,pml,2,0,1,0); - } - MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr); - tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period); - -#ifdef _OPENMP -#pragma omp parallel -#endif - for (MFIter mfi(reg); mfi.isValid(); ++mfi) + tmpregmf.setVal(0.0, 1, ncp-1, 0); + totpmlmf.ParallelCopy(tmpregmf,0, 0, 1, IntVect(0), ngp, period); +// #ifdef _OPENMP +// #pragma omp parallel +// #endif + for (MFIter mfi(pml); mfi.isValid(); ++mfi) { - const FArrayBox& src = tmpregmf[mfi]; - FArrayBox& dst = reg[mfi]; - const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); + const FArrayBox& src = totpmlmf[mfi]; + FArrayBox& dst = pml[mfi]; + const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); //amrex::boxDiff(dst.box(), mfi.validbox()); for (const Box& bx : bl) { dst.copy(src, bx, 0, bx, 0, 1); } } } +} + + +void +PML::CopyRegInPMLs (MultiFab& pml, MultiFab& reg, const Geometry& geom) +{ +// #ifdef _OPENMP +// #pragma omp parallel +// #endif + const IntVect& ngr = reg.nGrowVect(); + const IntVect& ngp = pml.nGrowVect(); + // const int ncp = pml.nComp(); + const auto& period = geom.periodicity(); + + // MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), 1, ngr); + // tmpregmf.setVal(0.0, 0, 1, ngr); + // MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, ngp); + // totpmlmf.setVal(0.0, 0, 1, ngp); + // realise sum of splitted fields inside pml + + + pml.ParallelCopy(reg, 0, 0, 1, ngr, ngp, period); + +} + +void +PML::CopyPMLsInReg (MultiFab& pml, MultiFab& reg, const Geometry& geom) +{ +// #ifdef _OPENMP +// #pragma omp parallel +// #endif + const IntVect& ngr = reg.nGrowVect(); + const IntVect& ngp = pml.nGrowVect(); + const auto& period = geom.periodicity(); + + // reg.ParallelCopy(pml, 0, 0, 1, IntVect(0), ngr, period); + reg.ParallelCopy(pml, 0, 0, 1, ngp, ngr, period); - // Copy from regular data to PML's first component - // Zero out the second (and third) component - MultiFab::Copy(tmpregmf,reg,0,0,1,0); - tmpregmf.setVal(0.0, 1, ncp-1, 0); - pml.ParallelCopy(tmpregmf, 0, 0, ncp, IntVect(0), ngp, period); } void diff --git a/Source/BoundaryConditions/PML_Cplx.cpp b/Source/BoundaryConditions/PML_Cplx.cpp new file mode 100644 index 000000000..0a47c5d84 --- /dev/null +++ b/Source/BoundaryConditions/PML_Cplx.cpp @@ -0,0 +1,958 @@ + +#include <PML.H> +#include <WarpX.H> +#include <WarpXConst.H> + +#include <AMReX_Print.H> +#include <AMReX_VisMF.H> + +#include <algorithm> + +#ifdef _OPENMP +#include <omp.h> +#endif + +using namespace amrex; + +namespace +{ + static void FillLo (int idim, Sigma& sigma, Sigma& sigma_star, + const Box& overlap, const Box& grid, Real fac) + { + int glo = grid.smallEnd(idim); + int olo = overlap.smallEnd(idim); + int ohi = overlap.bigEnd(idim); + int slo = sigma.m_lo; + int sslo = sigma_star.m_lo; + for (int i = olo; i <= ohi+1; ++i) + { + Real offset = static_cast<Real>(glo-i); + sigma[i-slo] = fac*(offset*offset); + } + for (int i = olo; i <= ohi; ++i) + { + Real offset = static_cast<Real>(glo-i) - 0.5; + sigma_star[i-sslo] = fac*(offset*offset); + } + } + + static void FillHi (int idim, Sigma& sigma, Sigma& sigma_star, + const Box& overlap, const Box& grid, Real fac) + { + int ghi = grid.bigEnd(idim); + int olo = overlap.smallEnd(idim); + int ohi = overlap.bigEnd(idim); + int slo = sigma.m_lo; + int sslo = sigma_star.m_lo; + for (int i = olo; i <= ohi+1; ++i) + { + Real offset = static_cast<Real>(i-ghi-1); + sigma[i-slo] = fac*(offset*offset); + } + for (int i = olo; i <= ohi; ++i) + { + Real offset = static_cast<Real>(i-ghi) - 0.5; + sigma_star[i-sslo] = fac*(offset*offset); + } + } + + static void FillZero (int idim, Sigma& sigma, Sigma& sigma_star, const Box& overlap) + { + int olo = overlap.smallEnd(idim); + int ohi = overlap.bigEnd(idim); + int slo = sigma.m_lo; + int sslo = sigma_star.m_lo; + std::fill(sigma.begin()+(olo-slo), sigma.begin()+(ohi+2-slo), 0.0); + std::fill(sigma_star.begin()+(olo-sslo), sigma_star.begin()+(ohi+1-sslo), 0.0); + } +} + +SigmaBox::SigmaBox (const Box& box, const BoxArray& grids, const Real* dx, int ncell, int delta) +{ + BL_ASSERT(box.cellCentered()); + + const IntVect& sz = box.size(); + const int* lo = box.loVect(); + const int* hi = box.hiVect(); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + sigma [idim].resize(sz[idim]+1); + sigma_star [idim].resize(sz[idim] ); + sigma_fac [idim].resize(sz[idim]+1); + sigma_star_fac[idim].resize(sz[idim] ); + + sigma [idim].m_lo = lo[idim]; + sigma [idim].m_hi = hi[idim]+1; + sigma_star [idim].m_lo = lo[idim]; + sigma_star [idim].m_hi = hi[idim]; + sigma_fac [idim].m_lo = lo[idim]; + sigma_fac [idim].m_hi = hi[idim]+1; + sigma_star_fac[idim].m_lo = lo[idim]; + sigma_star_fac[idim].m_hi = hi[idim]; + } + + Array<Real,AMREX_SPACEDIM> fac; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + fac[idim] = 4.0*PhysConst::c/(dx[idim]*static_cast<Real>(delta*delta)); + } + + const std::vector<std::pair<int,Box> >& isects = grids.intersections(box, false, ncell); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + int jdim = (idim+1) % AMREX_SPACEDIM; +#if (AMREX_SPACEDIM == 3) + int kdim = (idim+2) % AMREX_SPACEDIM; +#endif + + Vector<int> direct_faces, side_faces, direct_side_edges, side_side_edges, corners; + for (const auto& kv : isects) + { + const Box& grid_box = grids[kv.first]; + + if (amrex::grow(grid_box, idim, ncell).intersects(box)) + { + direct_faces.push_back(kv.first); + } + else if (amrex::grow(grid_box, jdim, ncell).intersects(box)) + { + side_faces.push_back(kv.first); + } +#if (AMREX_SPACEDIM == 3) + else if (amrex::grow(grid_box, kdim, ncell).intersects(box)) + { + side_faces.push_back(kv.first); + } + else if (amrex::grow(amrex::grow(grid_box,idim,ncell), + jdim,ncell).intersects(box)) + { + direct_side_edges.push_back(kv.first); + } + else if (amrex::grow(amrex::grow(grid_box,idim,ncell), + kdim,ncell).intersects(box)) + { + direct_side_edges.push_back(kv.first); + } + else if (amrex::grow(amrex::grow(grid_box,jdim,ncell), + kdim,ncell).intersects(box)) + { + side_side_edges.push_back(kv.first); + } +#endif + else + { + corners.push_back(kv.first); + } + } + + for (auto gid : corners) + { + const Box& grid_box = grids[gid]; + + Box lobox = amrex::adjCellLo(grid_box, idim, ncell); + lobox.grow(jdim,ncell); +#if (AMREX_SPACEDIM == 3) + lobox.grow(kdim,ncell); +#endif + Box looverlap = lobox & box; + if (looverlap.ok()) { + FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + } + + Box hibox = amrex::adjCellHi(grid_box, idim, ncell); + hibox.grow(jdim,ncell); +#if (AMREX_SPACEDIM == 3) + hibox.grow(kdim,ncell); +#endif + Box hioverlap = hibox & box; + if (hioverlap.ok()) { + FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + } + + if (!looverlap.ok() && !hioverlap.ok()) { + amrex::Abort("SigmaBox::SigmaBox(): corners, how did this happen?\n"); + } + } + +#if (AMREX_SPACEDIM == 3) + for (auto gid : side_side_edges) + { + const Box& grid_box = grids[gid]; + const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; + if (overlap.ok()) { + FillZero(idim, sigma[idim], sigma_star[idim], overlap); + } else { + amrex::Abort("SigmaBox::SigmaBox(): side_side_edges, how did this happen?\n"); + } + } + + for (auto gid : direct_side_edges) + { + const Box& grid_box = grids[gid]; + + Box lobox = amrex::adjCellLo(grid_box, idim, ncell); + Box looverlap = lobox.grow(jdim,ncell).grow(kdim,ncell) & box; + if (looverlap.ok()) { + FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + } + + Box hibox = amrex::adjCellHi(grid_box, idim, ncell); + Box hioverlap = hibox.grow(jdim,ncell).grow(kdim,ncell) & box; + if (hioverlap.ok()) { + FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + } + + if (!looverlap.ok() && !hioverlap.ok()) { + amrex::Abort("SigmaBox::SigmaBox(): direct_side_edges, how did this happen?\n"); + } + } +#endif + + for (auto gid : side_faces) + { + const Box& grid_box = grids[gid]; +#if (AMREX_SPACEDIM == 2) + const Box& overlap = amrex::grow(grid_box,jdim,ncell) & box; +#else + const Box& overlap = amrex::grow(amrex::grow(grid_box,jdim,ncell),kdim,ncell) & box; +#endif + if (overlap.ok()) { + FillZero(idim, sigma[idim], sigma_star[idim], overlap); + } else { + amrex::Abort("SigmaBox::SigmaBox(): side_faces, how did this happen?\n"); + } + } + + for (auto gid : direct_faces) + { + const Box& grid_box = grids[gid]; + + const Box& lobox = amrex::adjCellLo(grid_box, idim, ncell); + Box looverlap = lobox & box; + if (looverlap.ok()) { + FillLo(idim, sigma[idim], sigma_star[idim], looverlap, grid_box, fac[idim]); + } + + const Box& hibox = amrex::adjCellHi(grid_box, idim, ncell); + Box hioverlap = hibox & box; + if (hioverlap.ok()) { + FillHi(idim, sigma[idim], sigma_star[idim], hioverlap, grid_box, fac[idim]); + } + + if (!looverlap.ok() && !hioverlap.ok()) { + amrex::Abort("SigmaBox::SigmaBox(): direct faces, how did this happen?\n"); + } + } + + if (direct_faces.size() > 1) { + amrex::Abort("SigmaBox::SigmaBox(): direct_faces.size() > 1, Box gaps not wide enough?\n"); + } + } +} + +void +SigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) +{ + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + 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); + } + } + } +} + +void +SigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) +{ + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + 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); + } + } + } +} + +MultiSigmaBox::MultiSigmaBox (const BoxArray& ba, const DistributionMapping& dm, + const BoxArray& grid_ba, const Real* dx, int ncell, int delta) + : FabArray<SigmaBox>(ba,dm,1,0,MFInfo(), + FabFactory<SigmaBox>(grid_ba,dx,ncell,delta)) +{} + +void +MultiSigmaBox::ComputePMLFactorsB (const Real* dx, Real dt) +{ + if (dt == dt_B) return; + + dt_B = dt; + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(*this); mfi.isValid(); ++mfi) + { + (*this)[mfi].ComputePMLFactorsB(dx, dt); + } +} + +void +MultiSigmaBox::ComputePMLFactorsE (const Real* dx, Real dt) +{ + if (dt == dt_E) return; + + dt_E = dt; + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(*this); mfi.isValid(); ++mfi) + { + (*this)[mfi].ComputePMLFactorsE(dx, 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) + : m_geom(geom), + m_cgeom(cgeom) +{ + const BoxArray& ba = MakeBoxArray(*geom, grid_ba, ncell); + if (ba.size() == 0) { + m_ok = false; + return; + } else { + m_ok = true; + } + + 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 ); + + 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)); + pml_E_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Ez_nodal_flag), dm, 3, nge)); + pml_B_fp[0].reset(new MultiFab(amrex::convert(ba,WarpX::Bx_nodal_flag), dm, 2, ngb)); + pml_B_fp[1].reset(new MultiFab(amrex::convert(ba,WarpX::By_nodal_flag), dm, 2, ngb)); + pml_B_fp[2].reset(new MultiFab(amrex::convert(ba,WarpX::Bz_nodal_flag), dm, 2, ngb)); + + pml_E_fp[0]->setVal(0.0); + pml_E_fp[1]->setVal(0.0); + pml_E_fp[2]->setVal(0.0); + pml_B_fp[0]->setVal(0.0); + pml_B_fp[1]->setVal(0.0); + pml_B_fp[2]->setVal(0.0); + + if (do_dive_cleaning) + { + pml_F_fp.reset(new MultiFab(amrex::convert(ba,IntVect::TheUnitVector()), dm, 3, ngf)); + pml_F_fp->setVal(0.0); + } + + sigba_fp.reset(new MultiSigmaBox(ba, dm, grid_ba, geom->CellSize(), ncell, delta)); + + if (cgeom) + { + + nge = 1; + ngb = 1; + + BoxArray grid_cba = grid_ba; + grid_cba.coarsen(ref_ratio); + const BoxArray& cba = MakeBoxArray(*cgeom, grid_cba, ncell); + + DistributionMapping cdm{cba}; + + pml_E_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::Ex_nodal_flag), cdm, 3, nge)); + pml_E_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::Ey_nodal_flag), cdm, 3, nge)); + pml_E_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::Ez_nodal_flag), cdm, 3, nge)); + pml_B_cp[0].reset(new MultiFab(amrex::convert(cba,WarpX::Bx_nodal_flag), cdm, 2, ngb)); + pml_B_cp[1].reset(new MultiFab(amrex::convert(cba,WarpX::By_nodal_flag), cdm, 2, ngb)); + pml_B_cp[2].reset(new MultiFab(amrex::convert(cba,WarpX::Bz_nodal_flag), cdm, 2, ngb)); + + pml_E_cp[0]->setVal(0.0); + pml_E_cp[1]->setVal(0.0); + pml_E_cp[2]->setVal(0.0); + pml_B_cp[0]->setVal(0.0); + pml_B_cp[1]->setVal(0.0); + pml_B_cp[2]->setVal(0.0); + + if (do_dive_cleaning) + { + pml_F_cp.reset(new MultiFab(amrex::convert(cba,IntVect::TheUnitVector()), cdm, 3, ngf)); + pml_F_cp->setVal(0.0); + } + + sigba_cp.reset(new MultiSigmaBox(cba, cdm, grid_cba, cgeom->CellSize(), ncell, delta)); + } + +} + +// BoxArray +// PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell) +// { +// Box domain = geom.Domain(); +// for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { +// if ( ! geom.isPeriodic(idim) ) { +// domain.grow(idim, ncell); +// } +// } +// +// BoxList bl; +// for (int i = 0, N = grid_ba.size(); i < N; ++i) +// { +// const Box& grid_bx = grid_ba[i]; +// const IntVect& grid_bx_sz = grid_bx.size(); +// AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, +// "Consider using larger amr.blocking_factor"); +// +// Box bx = grid_bx; +// bx.grow(ncell); +// bx &= domain; +// +// Vector<Box> bndryboxes; +// #if (AMREX_SPACEDIM == 3) +// int kbegin = -1, kend = 1; +// #else +// int kbegin = 0, kend = 0; +// #endif +// for (int kk = kbegin; kk <= kend; ++kk) { +// for (int jj = -1; jj <= 1; ++jj) { +// for (int ii = -1; ii <= 1; ++ii) { +// if (ii != 0 || jj != 0 || kk != 0) { +// Box b = grid_bx; +// b.shift(grid_bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); +// b &= bx; +// if (b.ok()) { +// bndryboxes.push_back(b); +// } +// } +// } +// } +// } +// +// const BoxList& noncovered = grid_ba.complementIn(bx); +// for (const Box& b : noncovered) { +// for (const auto& bb : bndryboxes) { +// Box ib = b & bb; +// if (ib.ok()) { +// bl.push_back(ib); +// } +// } +// } +// } +// +// BoxArray ba(bl); +// ba.removeOverlap(false); +// +// return ba; +// } + +BoxArray +PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, int ncell) +{ + Box domain = geom.Domain(); + + IntVect limInfDomain = domain.smallEnd(); + IntVect limSupDomain = domain.bigEnd(); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( ! geom.isPeriodic(idim) ) { + // Create the simulation area (inside the whole domain) (only if non-periodical boundary conditions) + domain.grow(idim, -ncell); // I don't know if that works + // domain.growHi(idim, -ncell); + // domain.growLo(idim, -ncell); + } + } + + BoxList bl; + for (int i = 0, N = grid_ba.size(); i < N; ++i) + { + const Box& grid_bx = grid_ba[i]; + // const IntVect& grid_bx_sz = grid_bx.size(); + // AMREX_ALWAYS_ASSERT_WITH_MESSAGE(grid_bx.shortside() > ncell, + // "Consider using larger amr.blocking_factor"); + + Box bx = grid_bx; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + bx.grow(idim, -ncell); + } + bx &= domain; // is this step necessary? We're always inside the domain by doing so, even if there are some periodical boundary conditions... + const IntVect& bx_sz = bx.size(); + + Vector<Box> bndryboxes; +#if (AMREX_SPACEDIM == 3) + int kbegin = -1, kend = 1; +#else + int kbegin = 0, kend = 0; +#endif + for (int kk = kbegin; kk <= kend; ++kk) { + for (int jj = -1; jj <= 1; ++jj) { + for (int ii = -1; ii <= 1; ++ii) { + if (ii != 0 || jj != 0 || kk != 0) { + // if ((ii != 0 && jj == 0 && kk == 0)||(ii == 0 && jj != 0 && kk == 0)||(ii == 0 && jj == 0 && kk != 0)) { + if(ii!=0 && jj!=0){ + int xlimg, xlimd, ylimg, ylimd; + if (ii == -1){ + xlimg = grid_bx.smallEnd()[0]; + xlimd = limInfDomain[0]; + } + else if (ii = 1) { + xlimg = grid_bx.bigEnd()[0]; + xlimd = limSupDomain[0]; + } + + if (jj == -1){ + ylimg = grid_bx.smallEnd()[1]; + ylimd = limInfDomain[1]; + } + else if (jj = 1){ + ylimg = grid_bx.bigEnd()[1]; + ylimd = limSupDomain[1]; + } + + if (xlimd==xlimg && ylimd==ylimg){ + amrex::Print() << "INDICES : i = " << ii << " j = " << jj << std::endl; + Box b = grid_bx; //grid_bx; + b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b.shift(ncell * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b &= grid_bx;// grid_bx; + amrex::Print() << "BOITE = [" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "<<b.bigEnd()[0] << ", "<< b.bigEnd()[1] << "]," << std::endl; + if (b.ok()) { + bndryboxes.push_back(b); + amrex::Print() << "WAS SUCCESSFULY ADDED!" << std::endl; + } + } + + } + + // + // Box b; + // if ((xlimg==xlimd)||(ylimg==ylimd)){ + // b = grid_bx; //grid_bx; grid_bx = GRANDE BOITE + // b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + // b.shift(ncell * IntVect{AMREX_D_DECL(ii,jj,kk)}); + // } + // else { + // b = bx; //grid_bx; bx = PETITE BOITE + // b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + // } + + // b.shift(grid_bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + else if (ii == 0){ + int ylimg, ylimd; + if (jj == -1){ + ylimg = grid_bx.smallEnd()[1]; + ylimd = limInfDomain[1]; + } + else if (jj = 1){ + ylimg = grid_bx.bigEnd()[1]; + ylimd = limSupDomain[1]; + } + + if (ylimd==ylimg){ + Box b = grid_bx; //grid_bx; + b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b.shift(ncell * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b &= grid_bx;// grid_bx; + if (b.ok()) { + bndryboxes.push_back(b); + } + + } + + else { + Box b = bx; //grid_bx; + b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b &= grid_bx;// grid_bx; + if (b.ok()) { + bndryboxes.push_back(b); + } + } + + } + + else if (jj == 0){ + int xlimg, xlimd; + if (ii == -1){ + xlimg = grid_bx.smallEnd()[0]; + xlimd = limInfDomain[0]; + } + else if (ii = 1) { + xlimg = grid_bx.bigEnd()[0]; + xlimd = limSupDomain[0]; + } + if (xlimd==xlimg){ + Box b = grid_bx; //grid_bx; + b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b.shift(ncell * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b &= grid_bx;// grid_bx; + if (b.ok()) { + bndryboxes.push_back(b); + } + + } + + else { + Box b = bx; //grid_bx; + b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + b &= grid_bx;// grid_bx; + if (b.ok()) { + bndryboxes.push_back(b); + } + } + + } + // Box b = grid_bx; //grid_bx; + // b.shift(bx_sz * IntVect{AMREX_D_DECL(ii,jj,kk)}); + // b.shift(ncell * IntVect{AMREX_D_DECL(ii,jj,kk)}); + // b &= grid_bx;// grid_bx; + // if (b.ok()) { + // bndryboxes.push_back(b); + // } + } + + } + } + } + + // const BoxList& grid_ba_reduced = grid_ba.complementIn(bx); + BoxArray domain_ba = BoxArray(domain); + + const BoxList& noncovered = domain_ba.complementIn(grid_bx); //grid_ba.complementIn(bx); + for (const Box& b : noncovered) { + for (const auto& bb : bndryboxes) { + Box ib = b & bb; + if (ib.ok()) { + bl.push_back(ib); + } + } + } + } + amrex::Print() << "Printing PML boxes BEFORE cleaning" << std::endl; + amrex::Print() << "[" << std::endl; + for (const Box& b: bl) { + amrex::Print() << "[" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "<<b.bigEnd()[0] << ", "<< b.bigEnd()[1] << "]," << std::endl; + } + amrex::Print()<< "];" << std::endl; + + BoxArray ba(bl); + ba.removeOverlap(false); + + BoxList bl_2 = BoxList(ba); + + amrex::Print() << "Printing PML boxes AFTER cleaning" << std::endl; + amrex::Print() << "[" << std::endl; + for (const Box& b: bl_2) { + amrex::Print() << "[" << b.smallEnd()[0]<<", "<< b.smallEnd()[1]<< ", "<<b.bigEnd()[0] << ", "<< b.bigEnd()[1] << "]," << std::endl; + } + amrex::Print()<< "];" << std::endl; + + return ba; +} + +void +PML::ComputePMLFactors (amrex::Real dt) +{ + if (sigba_fp) { + sigba_fp->ComputePMLFactorsB(m_geom->CellSize(), dt); + sigba_fp->ComputePMLFactorsE(m_geom->CellSize(), dt); + } + if (sigba_cp) { + sigba_cp->ComputePMLFactorsB(m_cgeom->CellSize(), dt); + sigba_cp->ComputePMLFactorsE(m_cgeom->CellSize(), dt); + } +} + +std::array<MultiFab*,3> +PML::GetE_fp () +{ + return {pml_E_fp[0].get(), pml_E_fp[1].get(), pml_E_fp[2].get()}; +} + +std::array<MultiFab*,3> +PML::GetB_fp () +{ + return {pml_B_fp[0].get(), pml_B_fp[1].get(), pml_B_fp[2].get()}; +} + +std::array<MultiFab*,3> +PML::GetE_cp () +{ + return {pml_E_cp[0].get(), pml_E_cp[1].get(), pml_E_cp[2].get()}; +} + +std::array<MultiFab*,3> +PML::GetB_cp () +{ + return {pml_B_cp[0].get(), pml_B_cp[1].get(), pml_B_cp[2].get()}; +} + +MultiFab* +PML::GetF_fp () +{ + return pml_F_fp.get(); +} + +MultiFab* +PML::GetF_cp () +{ + return pml_F_cp.get(); +} + +void +PML::ExchangeB (const std::array<amrex::MultiFab*,3>& B_fp, + const std::array<amrex::MultiFab*,3>& B_cp) +{ + 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], *Bp[0], *m_geom); + Exchange(*pml_B_fp[1], *Bp[1], *m_geom); + Exchange(*pml_B_fp[2], *Bp[2], *m_geom); + } + else if (patch_type == PatchType::coarse && pml_B_cp[0] && Bp[0]) + { + 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); + } +} + +void +PML::ExchangeE (const std::array<amrex::MultiFab*,3>& E_fp, + const std::array<amrex::MultiFab*,3>& E_cp) +{ + ExchangeE(PatchType::fine, E_fp); + ExchangeE(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], *Ep[0], *m_geom); + Exchange(*pml_E_fp[1], *Ep[1], *m_geom); + Exchange(*pml_E_fp[2], *Ep[2], *m_geom); + } + else if (patch_type == PatchType::coarse && pml_E_cp[0] && Ep[0]) + { + 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) +{ + 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 +PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom) +{ + const IntVect& ngr = reg.nGrowVect(); + const IntVect& ngp = pml.nGrowVect(); + const int ncp = pml.nComp(); + const auto& period = geom.periodicity(); + + MultiFab tmpregmf(reg.boxArray(), reg.DistributionMap(), ncp, ngr); + + if (ngp.max() > 0) // Copy from pml to the ghost cells of regular data + { + MultiFab totpmlmf(pml.boxArray(), pml.DistributionMap(), 1, 0); + MultiFab::LinComb(totpmlmf, 1.0, pml, 0, 1.0, pml, 1, 0, 1, 0); + if (ncp == 3) { + MultiFab::Add(totpmlmf,pml,2,0,1,0); + } + + MultiFab::Copy(tmpregmf, reg, 0, 0, 1, ngr); + tmpregmf.ParallelCopy(totpmlmf, 0, 0, 1, IntVect(0), ngr, period); + +#ifdef _OPENMP +#pragma omp parallel +#endif + for (MFIter mfi(reg); mfi.isValid(); ++mfi) + { + const FArrayBox& src = tmpregmf[mfi]; + FArrayBox& dst = reg[mfi]; + const BoxList& bl = amrex::boxDiff(dst.box(), mfi.validbox()); + for (const Box& bx : bl) + { + dst.copy(src, bx, 0, bx, 0, 1); + } + } + } + + // Copy from regular data to PML's first component + // Zero out the second (and third) component + MultiFab::Copy(tmpregmf,reg,0,0,1,0); + tmpregmf.setVal(0.0, 1, ncp-1, 0); + pml.ParallelCopy(tmpregmf, 0, 0, ncp, IntVect(0), ngp, period); +} + +void +PML::FillBoundary () +{ + FillBoundaryE(); + FillBoundaryB(); + FillBoundaryF(); +} + +void +PML::FillBoundaryE () +{ + 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(); + Vector<MultiFab*> mf{pml_E_fp[0].get(),pml_E_fp[1].get(),pml_E_fp[2].get()}; + amrex::FillBoundary(mf, period); + } + else if (patch_type == PatchType::coarse && pml_E_cp[0] && pml_E_cp[0]->nGrowVect().max() > 0) + { + const auto& period = m_cgeom->periodicity(); + Vector<MultiFab*> mf{pml_E_cp[0].get(),pml_E_cp[1].get(),pml_E_cp[2].get()}; + amrex::FillBoundary(mf, period); + } +} + +void +PML::FillBoundaryB () +{ + 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(); + Vector<MultiFab*> mf{pml_B_fp[0].get(),pml_B_fp[1].get(),pml_B_fp[2].get()}; + amrex::FillBoundary(mf, period); + } + else if (patch_type == PatchType::coarse && pml_B_cp[0]) + { + const auto& period = m_cgeom->periodicity(); + Vector<MultiFab*> mf{pml_B_cp[0].get(),pml_B_cp[1].get(),pml_B_cp[2].get()}; + amrex::FillBoundary(mf, period); + } +} + +void +PML::FillBoundaryF () +{ + 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); + } + 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); + } +} + +void +PML::CheckPoint (const std::string& dir) const +{ + if (pml_E_fp[0]) + { + VisMF::Write(*pml_E_fp[0], dir+"_Ex_fp"); + VisMF::Write(*pml_E_fp[1], dir+"_Ey_fp"); + VisMF::Write(*pml_E_fp[2], dir+"_Ez_fp"); + VisMF::Write(*pml_B_fp[0], dir+"_Bx_fp"); + VisMF::Write(*pml_B_fp[1], dir+"_By_fp"); + VisMF::Write(*pml_B_fp[2], dir+"_Bz_fp"); + } + + if (pml_E_cp[0]) + { + VisMF::Write(*pml_E_cp[0], dir+"_Ex_cp"); + VisMF::Write(*pml_E_cp[1], dir+"_Ey_cp"); + VisMF::Write(*pml_E_cp[2], dir+"_Ez_cp"); + VisMF::Write(*pml_B_cp[0], dir+"_Bx_cp"); + VisMF::Write(*pml_B_cp[1], dir+"_By_cp"); + VisMF::Write(*pml_B_cp[2], dir+"_Bz_cp"); + } +} + +void +PML::Restart (const std::string& dir) +{ + if (pml_E_fp[0]) + { + VisMF::Read(*pml_E_fp[0], dir+"_Ex_fp"); + VisMF::Read(*pml_E_fp[1], dir+"_Ey_fp"); + VisMF::Read(*pml_E_fp[2], dir+"_Ez_fp"); + VisMF::Read(*pml_B_fp[0], dir+"_Bx_fp"); + VisMF::Read(*pml_B_fp[1], dir+"_By_fp"); + VisMF::Read(*pml_B_fp[2], dir+"_Bz_fp"); + } + + if (pml_E_cp[0]) + { + VisMF::Read(*pml_E_cp[0], dir+"_Ex_cp"); + VisMF::Read(*pml_E_cp[1], dir+"_Ey_cp"); + VisMF::Read(*pml_E_cp[2], dir+"_Ez_cp"); + VisMF::Read(*pml_B_cp[0], dir+"_Bx_cp"); + VisMF::Read(*pml_B_cp[1], dir+"_By_cp"); + VisMF::Read(*pml_B_cp[2], dir+"_Bz_cp"); + } +} diff --git a/Source/BoundaryConditions/PML_routines.F90 b/Source/BoundaryConditions/PML_routines.F90 index 380e52934..faf32dfc7 100644 --- a/Source/BoundaryConditions/PML_routines.F90 +++ b/Source/BoundaryConditions/PML_routines.F90 @@ -255,53 +255,295 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & + & jx, jxlo, jxhi, & + & jy, jylo, jyhi, & + & jz, jzlo, jzhi, & + & flag, & + & pml_type, & + & sigjx, sjxlo, sjxhi, & + & sigjy, sjylo, sjyhi, & + & sigjz, sjzlo, sjzhi, & + & sigsjx, ssjxlo, ssjxhi, & + & sigsjy, ssjylo, ssjyhi, & + & sigsjz, ssjzlo, ssjzhi, & + & mudt, & & dtsdx, dtsdy, dtsdz) & bind(c,name='warpx_push_pml_evec_3d') integer, intent(in) :: xlo(3), xhi(3), ylo(3), yhi(3), zlo(3), zhi(3), & Exlo(3), Exhi(3), Eylo(3), Eyhi(3), Ezlo(3), Ezhi(3), & - Bxlo(3), Bxhi(3), Bylo(3), Byhi(3), Bzlo(3), Bzhi(3) - real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz + Bxlo(3), Bxhi(3), Bylo(3), Byhi(3), Bzlo(3), Bzhi(3), & + jxlo(3), jxhi(3), jylo(3), jyhi(3), jzlo(3), jzhi(3), & + sjxlo, sjxhi, sjylo, sjyhi, sjzlo, sjzhi, & + ssjxlo, ssjxhi, ssjylo, ssjyhi, ssjzlo, ssjzhi, & + flag, pml_type(5) + real(amrex_real), intent(in ) :: mudt, dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),Exlo(3):Exhi(3),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),Eylo(3):Eyhi(3),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),Ezlo(3):Ezhi(3),2) real(amrex_real), intent(in ) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),Bxlo(3):Bxhi(3),2) real(amrex_real), intent(in ) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),Bylo(3):Byhi(3),2) real(amrex_real), intent(in ) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),Bzlo(3):Bzhi(3),2) + real(amrex_real), intent(in ) :: jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2),jxlo(3):jxhi(3)) + real(amrex_real), intent(in ) :: jy (jylo(1):jyhi(1),jylo(2):jyhi(2),jylo(3):jyhi(3)) + real(amrex_real), intent(in ) :: jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2),jzlo(3):jzhi(3)) + real(amrex_real), intent(in ) :: sigjx (sjxlo:sjxhi) + real(amrex_real), intent(in ) :: sigjy (sjylo:sjyhi) + real(amrex_real), intent(in ) :: sigjz (sjzlo:sjzhi) + real(amrex_real), intent(in ) :: sigsjx (ssjxlo:ssjxhi) + real(amrex_real), intent(in ) :: sigsjy (ssjylo:ssjyhi) + real(amrex_real), intent(in ) :: sigsjz (ssjzlo:ssjzhi) integer :: i, j, k + real(amrex_real) :: alpha_xy, alpha_xz, alpha_yx, alpha_yz, alpha_zx, alpha_zy + + if (flag == 0) then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & + & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2)) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & + & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2)) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i ,j ,k,1)+By(i ,j ,k,2) & + & -By(i-1,j ,k,1)-By(i-1,j ,k,2)) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) + end do + end do + end do + + else !flag = 1 + !!!!! DIRECT FACE + if (pml_type(1)==0) then + if (pml_type(2)==0) then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + alpha_xy = 0.5 + alpha_xz = 0.5 + Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & + & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2))& + & -mudt*alpha_xy*jx(i,j,k) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& + & -mudt*alpha_xz*jx(i,j,k) + end do + end do + end do - do k = xlo(3), xhi(3) - do j = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & - & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2)) - Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & - & -By(i,j ,k-1,1)-By(i,j ,k-1,2)) - end do - end do - end do + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + alpha_yx = 1. + alpha_yz = 0. + Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & + & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2))& + & -mudt*alpha_yx*jy(i,j,k) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& + & -mudt*alpha_yz*jy(i,j,k) + end do + end do + end do - do k = ylo(3), yhi(3) - do j = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & - & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2)) - Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & - & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2)) - end do - end do - end do + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + alpha_zx = 1. + alpha_zy = 0. + Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i ,j ,k,1)+By(i ,j ,k,2) & + & -By(i-1,j ,k,1)-By(i-1,j ,k,2))& + & -mudt*alpha_zx*jz(i,j,k) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& + & -mudt*alpha_zy*jz(i,j,k) + end do + end do + end do - do k = zlo(3), zhi(3) - do j = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i ,j ,k,1)+By(i ,j ,k,2) & - & -By(i-1,j ,k,1)-By(i-1,j ,k,2)) - Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & - & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2)) - end do - end do - end do + else if (pml_type(2)==1) then + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + alpha_xy = 1. + alpha_xz = 0. + Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & + & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2))& + & -mudt*alpha_xy*jx(i,j,k) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& + & -mudt*alpha_xz*jx(i,j,k) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + alpha_yx = 0.5 + alpha_yz = 0.5 + Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & + & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2))& + & -mudt*alpha_yx*jy(i,j,k) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& + & -mudt*alpha_yz*jy(i,j,k) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + alpha_zx = 0. + alpha_zy = 1. + Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i ,j ,k,1)+By(i ,j ,k,2) & + & -By(i-1,j ,k,1)-By(i-1,j ,k,2))& + & -mudt*alpha_zx*jz(i,j,k) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& + & -mudt*alpha_zy*jz(i,j,k) + end do + end do + end do + + else !(pml_type(2)==2) + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + alpha_xy = 0. + alpha_xz = 1. + Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & + & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2))& + & -mudt*alpha_xy*jx(i,j,k) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& + & -mudt*alpha_xz*jx(i,j,k) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + alpha_yx = 0. + alpha_yz = 1. + Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & + & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2))& + & -mudt*alpha_yx*jy(i,j,k) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& + & -mudt*alpha_yz*jy(i,j,k) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + alpha_zx = 0.5 + alpha_zy = 0.5 + Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i ,j ,k,1)+By(i ,j ,k,2) & + & -By(i-1,j ,k,1)-By(i-1,j ,k,2))& + & -mudt*alpha_zx*jz(i,j,k) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& + & -mudt*alpha_zy*jz(i,j,k) + end do + end do + end do + end if !DIRECT FACE + + !!!!! SIDE EDGE OR CORNER + else + do k = xlo(3), xhi(3) + do j = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + ! compute current coefficients alpha + if ((sigjy(j)==0.) .AND. (sigjz(k)==0.)) then + alpha_xy = 0.5 + alpha_xz = 0.5 + else + alpha_xy = sigjy(j)/(sigjy(j)+sigjz(k)) + alpha_xz = sigjz(k)/(sigjy(j)+sigjz(k)) + end if + + Ex(i,j,k,1) = Ex(i,j,k,1) + dtsdy*(Bz(i,j ,k ,1)+Bz(i,j ,k ,2) & + & -Bz(i,j-1,k ,1)-Bz(i,j-1,k ,2))& + & -mudt*alpha_xy*jx(i,j,k) + Ex(i,j,k,2) = Ex(i,j,k,2) - dtsdz*(By(i,j ,k ,1)+By(i,j ,k ,2) & + & -By(i,j ,k-1,1)-By(i,j ,k-1,2))& + & -mudt*alpha_xz*jx(i,j,k) + end do + end do + end do + + do k = ylo(3), yhi(3) + do j = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + if ((sigjx(i)==0.) .AND. (sigjz(k)==0.)) then + alpha_yx = 0.5 + alpha_yz = 0.5 + else + alpha_yx = sigjx(i)/(sigjx(i)+sigjz(k)) + alpha_yz = sigjz(k)/(sigjx(i)+sigjz(k)) + end if + + Ey(i,j,k,1) = Ey(i,j,k,1) + dtsdz*(Bx(i ,j,k ,1)+Bx(i ,j,k ,2) & + & -Bx(i ,j,k-1,1)-Bx(i ,j,k-1,2))& + & -mudt*alpha_yx*jy(i,j,k) + Ey(i,j,k,2) = Ey(i,j,k,2) - dtsdx*(Bz(i ,j,k ,1)+Bz(i ,j,k ,2) & + & -Bz(i-1,j,k ,1)-Bz(i-1,j,k ,2))& + & -mudt*alpha_yz*jy(i,j,k) + end do + end do + end do + + do k = zlo(3), zhi(3) + do j = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + if ((sigjx(i)==0.) .AND. (sigjy(j)==0.)) then + alpha_zx = 0.5 + alpha_zy = 0.5 + else + alpha_zx = sigjx(i)/(sigjx(i)+sigjy(j)) + alpha_zy = sigjy(j)/(sigjx(i)+sigjy(j)) + end if + + Ez(i,j,k,1) = Ez(i,j,k,1) + dtsdx*(By(i ,j ,k,1)+By(i ,j ,k,2) & + & -By(i-1,j ,k,1)-By(i-1,j ,k,2))& + & -mudt*alpha_zx*jz(i,j,k) + Ez(i,j,k,2) = Ez(i,j,k,2) - dtsdy*(Bx(i ,j ,k,1)+Bx(i ,j ,k,2) & + & -Bx(i ,j-1,k,1)-Bx(i ,j-1,k,2))& + & -mudt*alpha_zy*jz(i,j,k) + end do + end do + end do + + end if !PML_TYPE + + end if ! FLAG end subroutine warpx_push_pml_evec_3d @@ -431,43 +673,83 @@ contains & Bx, Bxlo, Bxhi, & & By, Bylo, Byhi, & & Bz, Bzlo, Bzhi, & - & dtsdx, dtsdy, dtsdz) & + & jx, jxlo, jxhi, & + & jy, jylo, jyhi, & + & jz, jzlo, jzhi, & + & flag, & + & mudt, dtsdx, dtsdy, dtsdz) & bind(c,name='warpx_push_pml_evec_2d') integer, intent(in) :: xlo(2), xhi(2), ylo(2), yhi(2), zlo(2), zhi(2), & Exlo(2), Exhi(2), Eylo(2), Eyhi(2), Ezlo(2), Ezhi(2), & - Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2) - real(amrex_real), intent(in) :: dtsdx, dtsdy, dtsdz + Bxlo(2), Bxhi(2), Bylo(2), Byhi(2), Bzlo(2), Bzhi(2), & + jxlo(2), jxhi(2), jylo(2), jyhi(2), jzlo(2), jzhi(2), & + flag + real(amrex_real), intent(in) :: mudt, dtsdx, dtsdy, dtsdz real(amrex_real), intent(inout) :: Ex (Exlo(1):Exhi(1),Exlo(2):Exhi(2),2) real(amrex_real), intent(inout) :: Ey (Eylo(1):Eyhi(1),Eylo(2):Eyhi(2),2) real(amrex_real), intent(inout) :: Ez (Ezlo(1):Ezhi(1),Ezlo(2):Ezhi(2),2) real(amrex_real), intent(in ) :: Bx (Bxlo(1):Bxhi(1),Bxlo(2):Bxhi(2),2) real(amrex_real), intent(in ) :: By (Bylo(1):Byhi(1),Bylo(2):Byhi(2),2) real(amrex_real), intent(in ) :: Bz (Bzlo(1):Bzhi(1),Bzlo(2):Bzhi(2),2) + real(amrex_real), intent(in ) :: jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2)) !jx (jxlo(1):jxhi(1),jxlo(2):jxhi(2),1) + real(amrex_real), intent(in ) :: jy (jylo(1):jyhi(1),jylo(2):jyhi(2)) !jy (jylo(1):jyhi(1),jylo(2):jyhi(2),1) + real(amrex_real), intent(in ) :: jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2)) !jz (jzlo(1):jzhi(1),jzlo(2):jzhi(2),1) integer :: i, k - do k = xlo(2), xhi(2) - do i = xlo(1), xhi(1) - Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & - & -By(i,k-1,1)-By(i,k-1,2)) - end do - end do - - do k = ylo(2), yhi(2) - do i = ylo(1), yhi(1) - Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & - & -Bx(i ,k-1,1)-Bx(i ,k-1,2)) - Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & - & -Bz(i-1,k ,1)-Bz(i-1,k ,2)) - end do - end do - - do k = zlo(2), zhi(2) - do i = zlo(1), zhi(1) - Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & - & -By(i-1,k,1)-By(i-1,k,2)) - end do - end do + if (flag == 1) then + do k = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & + & -By(i,k-1,1)-By(i,k-1,2))& + & - mudt * jx(i,k) + end do + end do + + do k = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & + & -Bx(i ,k-1,1)-Bx(i ,k-1,2))& + & - mudt * 0.5 * jy(i,k) + Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & + & -Bz(i-1,k ,1)-Bz(i-1,k ,2))& + & - mudt * 0.5 * jy(i,k) + end do + end do + + do k = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & + & -By(i-1,k,1)-By(i-1,k,2))& + & - mudt * jz(i,k) + end do + end do + + else !no particles in PML + do k = xlo(2), xhi(2) + do i = xlo(1), xhi(1) + Ex(i,k,2) = Ex(i,k,2) - dtsdz*(By(i,k ,1)+By(i,k ,2) & + & -By(i,k-1,1)-By(i,k-1,2)) + end do + end do + + do k = ylo(2), yhi(2) + do i = ylo(1), yhi(1) + Ey(i,k,1) = Ey(i,k,1) + dtsdz*(Bx(i ,k ,1)+Bx(i ,k ,2) & + & -Bx(i ,k-1,1)-Bx(i ,k-1,2)) + Ey(i,k,2) = Ey(i,k,2) - dtsdx*(Bz(i ,k ,1)+Bz(i ,k ,2) & + & -Bz(i-1,k ,1)-Bz(i-1,k ,2)) + end do + end do + + do k = zlo(2), zhi(2) + do i = zlo(1), zhi(1) + Ez(i,k,1) = Ez(i,k,1) + dtsdx*(By(i ,k,1)+By(i ,k,2) & + & -By(i-1,k,1)-By(i-1,k,2)) + end do + end do + + end if !flag end subroutine warpx_push_pml_evec_2d @@ -827,6 +1109,7 @@ contains do k = texlo(2), texhi(2) do i = texlo(1), texhi(1) + ! ex(i,k,1) = ex(i,k,1) * sigez(k) !No damping for Exy ex(i,k,2) = ex(i,k,2) * sigez(k) ex(i,k,3) = ex(i,k,3) * sigcx(i) end do @@ -867,6 +1150,91 @@ contains end subroutine warpx_damp_pml_2d + subroutine warpx_dampJ_pml_2d (tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + & jx, jxlo, jxhi, jy, jylo, jyhi, jz, jzlo, jzhi, & + & sigjx, sjxlo, sjxhi, sigjz, sjzlo, sjzhi, & + & sigsjx, ssjxlo, ssjxhi, sigsjz, ssjzlo, ssjzhi) & + bind(c,name='warpx_dampJ_pml_2d') + integer, dimension(2), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + jxlo, jxhi, jylo, jyhi, jzlo, jzhi + integer, intent(in), value :: sjxlo, sjxhi, sjzlo, sjzhi, ssjxlo, ssjxhi, ssjzlo, ssjzhi + real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2)) + real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2)) + real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2)) + real(amrex_real), intent(in) :: sigjx(sjxlo:sjxhi) + real(amrex_real), intent(in) :: sigjz(sjzlo:sjzhi) + real(amrex_real), intent(in) :: sigsjx(ssjxlo:ssjxhi) + real(amrex_real), intent(in) :: sigsjz(ssjzlo:ssjzhi) + + integer :: i,k + + do k = tjxlo(2), tjxhi(2) + do i = tjxlo(1), tjxhi(1) + jx(i,k) = jx(i,k) * sigsjx(i) * sigjz(k) + end do + end do + + ! do k = tjylo(2), tjyhi(2) + ! do i = tjylo(1), tjyhi(1) + ! jy(i,k) = jy(i,k) + ! end do + ! end do + + do k = tjzlo(2), tjzhi(2) + do i = tjzlo(1), tjzhi(1) + jz(i,k) = jz(i,k) * sigjx(i) * sigsjz(k) + end do + end do + + + end subroutine warpx_dampJ_pml_2d + + subroutine warpx_dampJ_pml_3d (tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + & jx, jxlo, jxhi, jy, jylo, jyhi, jz, jzlo, jzhi, & + & sigjx, sjxlo, sjxhi, sigjy, sjylo, sjyhi, sigjz, sjzlo, sjzhi, & + & sigsjx, ssjxlo, ssjxhi, sigsjy, ssjylo, ssjyhi, sigsjz, ssjzlo, ssjzhi) & + bind(c,name='warpx_dampJ_pml_3d') + integer, dimension(3), intent(in) :: tjxlo, tjxhi, tjylo, tjyhi, tjzlo, tjzhi, & + jxlo, jxhi, jylo, jyhi, jzlo, jzhi + integer, intent(in), value :: sjxlo, sjxhi, sjylo, sjyhi, sjzlo, sjzhi, ssjxlo, ssjxhi, ssjylo, ssjyhi, ssjzlo, ssjzhi + real(amrex_real), intent(inout) :: jx(jxlo(1):jxhi(1),jxlo(2):jxhi(2),jxlo(3):jxhi(3)) + real(amrex_real), intent(inout) :: jy(jylo(1):jyhi(1),jylo(2):jyhi(2),jylo(3):jyhi(3)) + real(amrex_real), intent(inout) :: jz(jzlo(1):jzhi(1),jzlo(2):jzhi(2),jzlo(3):jzhi(3)) + real(amrex_real), intent(in) :: sigjx(sjxlo:sjxhi) + real(amrex_real), intent(in) :: sigjy(sjylo:sjyhi) + real(amrex_real), intent(in) :: sigjz(sjzlo:sjzhi) + real(amrex_real), intent(in) :: sigsjx(ssjxlo:ssjxhi) + real(amrex_real), intent(in) :: sigsjy(ssjylo:ssjyhi) + real(amrex_real), intent(in) :: sigsjz(ssjzlo:ssjzhi) + + integer :: i,j,k + + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jx(i,j,k) = jx(i,j,k) * sigsjx(i) * sigjy(j) * sigjz(k) + end do + end do + end do + + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jy(i,j,k) = jy(i,j,k) * sigjx(i) * sigsjy(j) * sigjz(k) + end do + end do + end do + + do k = tjxlo(3), tjxhi(3) + do j = tjylo(2), tjyhi(2) + do i = tjxlo(1), tjxhi(1) + jz(i,j,k) = jz(i,j,k) * sigjx(i) * sigjy(j) * sigsjz(k) + end do + end do + end do + + end subroutine warpx_dampJ_pml_3d + subroutine warpx_damp_pml_3d (texlo, texhi, teylo, teyhi, tezlo, tezhi, & & tbxlo, tbxhi, tbylo, tbyhi, tbzlo, tbzhi, & diff --git a/Source/BoundaryConditions/WarpXEvolvePML.cpp b/Source/BoundaryConditions/WarpXEvolvePML.cpp index b0688b2c1..08173e424 100644 --- a/Source/BoundaryConditions/WarpXEvolvePML.cpp +++ b/Source/BoundaryConditions/WarpXEvolvePML.cpp @@ -55,7 +55,6 @@ WarpX::DampPML (int lev, PatchType patch_type) 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(), @@ -79,3 +78,55 @@ WarpX::DampPML (int lev, PatchType patch_type) } } } + +void +WarpX::DampJPML () +{ + for (int lev = 0; lev <= finest_level; ++lev) { + DampJPML(lev); + } +} + +void +WarpX::DampJPML (int lev) +{ + DampJPML(lev, PatchType::fine); + if (lev > 0) DampJPML(lev, PatchType::coarse); +} + +void +WarpX::DampJPML (int lev, PatchType patch_type) +{ + if (!do_pml) return; + if (!do_pml_j_damping) return; + + BL_PROFILE("WarpX::DampJPML()"); + + if (pml[lev]->ok()) + { + const auto& pml_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); + const auto& sigba = (patch_type == PatchType::fine) ? pml[lev]->GetMultiSigmaBox_fp() + : pml[lev]->GetMultiSigmaBox_cp(); + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + // amrex::Print()<<"===== DAMP PML ====="<<std::endl; + for ( MFIter mfi(*pml_j[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + const Box& tjx = mfi.tilebox(jx_nodal_flag); + const Box& tjy = mfi.tilebox(jy_nodal_flag); + const Box& tjz = mfi.tilebox(jz_nodal_flag); + // amrex::Print()<< "tjx = ["<< tjx.smallEnd()[0]<<", "<<tjx.smallEnd()[1]<<", "<<tjx.bigEnd()[0]<<", "<<tjx.bigEnd()[1]<<"]" <<std::endl; + // amrex::Print()<< "tjz = ["<< tjz.smallEnd()[0]<<", "<<tjz.smallEnd()[1]<<", "<<tjz.bigEnd()[0]<<", "<<tjz.bigEnd()[1]<<"]" <<std::endl; + WRPX_DAMPJ_PML(tjx.loVect(), tjx.hiVect(), + tjy.loVect(), tjy.hiVect(), + tjz.loVect(), tjz.hiVect(), + BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), + WRPX_PML_SIGMAJ_TO_FORTRAN(sigba[mfi])); //WRPX_PML_SIGMACOEFF_TO_FORTRAN(sigba[mfi]) + } + + } +} diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 32a4747db..531cd6a56 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -135,7 +135,7 @@ WarpX::EvolveEM (int numsteps) bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); // slice generation // - bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0); + bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0); bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); @@ -145,7 +145,7 @@ WarpX::EvolveEM (int numsteps) // We might need to move j because we are going to make a plotfile. int num_moved = MoveWindow(move_j); - + if (max_level == 0) { int num_redistribute_ghost = num_moved + 1; mypc->RedistributeLocal(num_redistribute_ghost); @@ -228,7 +228,7 @@ WarpX::EvolveEM (int numsteps) // End loop on time steps } - bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step + bool write_plot_file = plot_int > 0 && istep[0] > last_plot_file_step && (max_time_reached || istep[0] >= max_step); bool do_insitu = (insitu_start >= istep[0]) && (insitu_int > 0) && @@ -255,7 +255,7 @@ WarpX::EvolveEM (int numsteps) UpdateInSitu(); } - if (check_int > 0 && istep[0] > last_check_file_step && + if (check_int > 0 && istep[0] > last_check_file_step && (max_time_reached || istep[0] >= max_step)) { WriteCheckPointFile(); } @@ -302,19 +302,75 @@ WarpX::OneStep_nosub (Real cur_time) FillBoundaryE(); FillBoundaryB(); #else + // amrex::Print()<< "WarpXEvolveEM.cpp : before CopyJinPMLs "<<std::endl; + if (do_pml && pml_has_particles){ + // do current deposition in PMLs + // copy current computed on mother grid to PMLs + // amrex::Print()<< "WarpXEvolveEM.cpp : IN CopyJinPMLs "<<std::endl; + for (int lev = 0; lev <= finest_level; ++lev) + { + if (pml[lev]->ok()){ + pml[lev]->CopyJinPMLs({ current_fp[lev][0].get(), + current_fp[lev][1].get(), + current_fp[lev][2].get() }, + { current_cp[lev][0].get(), + current_cp[lev][1].get(), + current_cp[lev][2].get() }); + } + } + } + + if (do_pml && do_pml_j_damping){ + // amrex::Print()<< "WarpXEvolveEM.cpp : DampJ "<<std::endl; + // damp current in pmls + // amrex::Print() << "===== DAMPING IN PMLs =====" << std::endl; + + // amrex::Print()<< "===== DAMPING J ====="<< std::endl; + // amrex::Print()<< "[AV DAMP] max_Jx_pml = "<< pml[0]->Getj_fp()[0]->min(0) << std::endl; + DampJPML(); + // amrex::Print()<< "[AP DAMP] max_Jx_pml = "<< pml[0]->Getj_fp()[0]->min(0) << std::endl; + } + EvolveF(0.5*dt[0], DtType::FirstHalf); FillBoundaryF(); EvolveB(0.5*dt[0]); // We now have B^{n+1/2} FillBoundaryB(); + EvolveE(dt[0]); // We now have E^{n+1} FillBoundaryE(); EvolveF(0.5*dt[0], DtType::SecondHalf); + + // amrex::Print()<< "WarpXEvolveEM.cpp : before CopyJinReg "<<std::endl; + if (do_pml && pml_has_particles){ + // do current deposition in PMLs + // copy current computed on mother grid to PMLs + // amrex::Print()<< "WarpXEvolveEM.cpp : IN CopyJinReg "<<std::endl; + for (int lev = 0; lev <= finest_level; ++lev) + { + if (pml[lev]->ok()){ + // amrex::Print()<< "[AV COPY] max_Jx = "<< current_fp[lev][0].get()->min(0) << std::endl; + // amrex::Print()<< "[AV COPY] max_Jx_pml = "<< pml[lev]->Getj_fp()[0]->min(0) << std::endl; + // amrex::Print()<< "===== Copy J from PML to Reg ====="<< std::endl; + pml[lev]->CopyJinReg({ current_fp[lev][0].get(), + current_fp[lev][1].get(), + current_fp[lev][2].get() }, + { current_cp[lev][0].get(), + current_cp[lev][1].get(), + current_cp[lev][2].get() }); + // amrex::Print()<< "[AP COPY] max_Jx = "<< current_fp[lev][0].get()->min(0) << std::endl; + + } + } + } + EvolveB(0.5*dt[0]); // We now have B^{n+1} if (do_pml) { + // amrex::Print()<< "WarpXEvolveEM.cpp : Damp "<<std::endl; DampPML(); FillBoundaryE(); } FillBoundaryB(); + #endif } @@ -524,13 +580,13 @@ WarpX::ComputeDt () /* \brief computes max_step for wakefield simulation in boosted frame. * \param geom: Geometry object that contains simulation domain. - * - * max_step is set so that the simulation stop when the lower corner of the + * + * max_step is set so that the simulation stop when the lower corner of the * simulation box passes input parameter zmax_plasma_to_compute_max_step. */ void WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ - // Sanity checks: can use zmax_plasma_to_compute_max_step only if + // Sanity checks: can use zmax_plasma_to_compute_max_step only if // the moving window and the boost are all in z direction. AMREX_ALWAYS_ASSERT_WITH_MESSAGE( WarpX::moving_window_dir == AMREX_SPACEDIM-1, @@ -547,7 +603,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ "Can use zmax_plasma_to_compute_max_step only if " + "warpx.boost_direction = z. TODO: all directions."); - // Lower end of the simulation domain. All quantities are given in boosted + // Lower end of the simulation domain. All quantities are given in boosted // frame except zmax_plasma_to_compute_max_step. const Real zmin_domain_boost = a_geom.ProbLo(AMREX_SPACEDIM-1); // End of the plasma: Transform input argument @@ -555,7 +611,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; // Plasma velocity const Real v_plasma_boost = -beta_boost * PhysConst::c; - // Get time at which the lower end of the simulation domain passes the + // Get time at which the lower end of the simulation domain passes the // upper end of the plasma (in the z direction). const Real interaction_time_boost = (len_plasma_boost-zmin_domain_boost)/ (moving_window_v-v_plasma_boost); @@ -568,7 +624,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ /* \brief Apply perfect mirror condition inside the box (not at a boundary). * In practice, set all fields to 0 on a section of the simulation domain - * (as for a perfect conductor with a given thickness). + * (as for a perfect conductor with a given thickness). * The mirror normal direction has to be parallel to the z axis. */ void @@ -585,10 +641,10 @@ WarpX::applyMirrors(Real time){ } // Loop over levels for(int lev=0; lev<=finest_level; lev++){ - // Make sure that the mirror contains at least + // Make sure that the mirror contains at least // mirror_z_npoints[i_mirror] cells Real dz = WarpX::CellSize(lev)[2]; - Real z_max = std::max(z_max_tmp, + Real z_max = std::max(z_max_tmp, z_min+mirror_z_npoints[i_mirror]*dz); // Get fine patch field MultiFabs MultiFab& Ex = *Efield_fp[lev][0].get(); diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index c53e13f8f..298b03dc6 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -137,6 +137,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { 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_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_cp(); #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -318,7 +319,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) 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_j = (patch_type == PatchType::fine) ? pml[lev]->Getj_fp() : pml[lev]->Getj_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 if (Gpu::notInLaunchRegion()) #endif @@ -328,6 +332,15 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) const Box& tey = mfi.tilebox(Ey_nodal_flag); const Box& tez = mfi.tilebox(Ez_nodal_flag); + // amrex::Print()<< "== sigba.pml_type_array"<<std::endl; + const auto& pml_type = sigba[mfi].pml_type_array; + // amrex::Print()<<"&pml_type = "<< &pml_type <<std::endl; + // amrex::Print()<<"pml_type[0] = "<< pml_type[0] <<std::endl; + // amrex::Print()<<"tex.loVect() = "<< tex.loVect() <<std::endl; + // amrex::Print()<<"sigba[mfi].sigma[0].m_lo = "<< sigba[mfi].sigma[0].m_lo <<std::endl; + // amrex::Print()<<"sigba[mfi].sigma[0].m_loVec() = "<< sigba[mfi].sigma[0].lo <<std::endl; + + // amrex::Print()<<"===== WRPX_PUSH_PML_EVEC ====="<<std::endl; WRPX_PUSH_PML_EVEC( tex.loVect(), tex.hiVect(), tey.loVect(), tey.hiVect(), @@ -338,8 +351,24 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) BL_TO_FORTRAN_3D((*pml_B[0])[mfi]), BL_TO_FORTRAN_3D((*pml_B[1])[mfi]), BL_TO_FORTRAN_3D((*pml_B[2])[mfi]), +#if (AMREX_SPACEDIM==2) + BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), + &pml_has_particles, + &mu_c2_dt, +#endif +#if (AMREX_SPACEDIM==3) + BL_TO_FORTRAN_3D((*pml_j[0])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[1])[mfi]), + BL_TO_FORTRAN_3D((*pml_j[2])[mfi]), + &pml_has_particles, + &pml_type, //sigba[mfi].pml_type_array, + WRPX_PML_SIGMACOEFF_TO_FORTRAN(sigba[mfi]), //WRPX_PML_SIGMAJ_TO_FORTRAN( + &mu_c2_dt, +#endif &dtsdx_c2, &dtsdy_c2, &dtsdz_c2); - + // amrex::Print()<<"===== WRPX_PUSH_PML_EVEC FIN ====="<<std::endl; if (pml_F) { WRPX_PUSH_PML_EVEC_F( @@ -435,4 +464,3 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) } } } - diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 98dd8685d..5049dce0e 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -38,6 +38,7 @@ #define WRPX_INTERPOLATE_CIC_TWO_LEVELS warpx_interpolate_cic_two_levels_3d #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_3d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_3d +#define WRPX_DAMPJ_PML warpx_dampJ_pml_3d #elif (AMREX_SPACEDIM == 2) @@ -51,6 +52,7 @@ #define WRPX_PUSH_PML_F warpx_push_pml_f_2d #define WRPX_DAMP_PML warpx_damp_pml_2d #define WRPX_DAMP_PML_F warpx_damp_pml_f_2d +#define WRPX_DAMPJ_PML warpx_dampJ_pml_2d #define WRPX_SUM_FINE_TO_CRSE_NODAL warpx_sum_fine_to_crse_nodal_2d #define WRPX_ZERO_OUT_BNDRY warpx_zero_out_bndry_2d @@ -377,6 +379,27 @@ extern "C" const BL_FORT_FAB_ARG_3D(bx), const BL_FORT_FAB_ARG_3D(by), const BL_FORT_FAB_ARG_3D(bz), +#if (AMREX_SPACEDIM == 2) + BL_FORT_FAB_ARG_3D(jx), + BL_FORT_FAB_ARG_3D(jy), + BL_FORT_FAB_ARG_3D(jz), + const int* flag, + const amrex::Real* mudt, +#endif +#if (AMREX_SPACEDIM == 3) + BL_FORT_FAB_ARG_3D(jx), + BL_FORT_FAB_ARG_3D(jy), + BL_FORT_FAB_ARG_3D(jz), + const int* flag, + const std::array<int,5>* pml_type,//const std::string pml_type, + const amrex::Real* sigjx, const int* sigjx_lo, const int* sigjx_hi, + const amrex::Real* sigjy, const int* sigjy_lo, const int* sigjy_hi, + const amrex::Real* sigjz, const int* sigjz_lo, const int* sigjz_hi, + const amrex::Real* sigsjx, const int* sigsjx_lo, const int* sigsjx_hi, + const amrex::Real* sigsjy, const int* sigsjy_lo, const int* sigsjy_hi, + const amrex::Real* sigjsz, const int* sigsjz_lo, const int* sigsjz_hi, + const amrex::Real* mudt, +#endif const amrex::Real* dtsdx, const amrex::Real* dtsdy, const amrex::Real* dtsdz); @@ -438,6 +461,25 @@ extern "C" #endif const amrex::Real* sigbz, int sigbz_lo, int sigbz_hi); + + void WRPX_DAMPJ_PML (const int* tjxlo, const int* tjxhi, + const int* tjylo, const int* tjyhi, + const int* tjzlo, const int* tjzhi, + amrex::Real* jx, const int* jxlo, const int* jxhi, + amrex::Real* jy, const int* jylo, const int* jyhi, + amrex::Real* jz, const int* jzlo, const int* jzhi, + const amrex::Real* sigjx, int sigjx_lo, int sigjx_hi, +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigjy, int sigjy_lo, int sigjy_hi, +#endif + const amrex::Real* sigjz, int sigjz_lo, int sigjz_hi, + const amrex::Real* sigsjx, int sigsjx_lo, int sigsjx_hi, +#if (AMREX_SPACEDIM == 3) + const amrex::Real* sigsjy, int sigsjy_lo, int sigsjy_hi, +#endif + const amrex::Real* sigjsz, int sigsjz_lo, int sigsjz_hi); + + void WRPX_SYNC_CURRENT (const int* lo, const int* hi, BL_FORT_FAB_ARG_ANYD(crse), const BL_FORT_FAB_ARG_ANYD(fine), diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index d583b2b0f..76704bcfe 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -88,7 +88,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 ); @@ -97,7 +97,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])); } } @@ -118,10 +118,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); } @@ -134,13 +134,13 @@ 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, do_dive_cleaning, do_moving_window, pml_has_particles)); //pml_has_particles)); 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)); + do_moving_window, pml_has_particles)); //pml_has_particles)); } } } @@ -226,7 +226,7 @@ WarpX::InitOpenbc () Vector<int> 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) { @@ -252,7 +252,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); diff --git a/Source/WarpX.H b/Source/WarpX.H index 4ad3d119f..0587e82d2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -149,12 +149,12 @@ public: BilinearFilter bilinear_filter; amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_exeybz; amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_bxbyez; - + static int num_mirrors; amrex::Vector<amrex::Real> mirror_z; amrex::Vector<amrex::Real> mirror_z_width; amrex::Vector<int> mirror_z_npoints; - + void applyMirrors(amrex::Real time); void ComputeDt (); @@ -179,6 +179,10 @@ public: void DampPML (int lev); void DampPML (int lev, PatchType patch_type); + void DampJPML (); + void DampJPML (int lev); + void DampJPML (int lev, PatchType patch_type); + void PushParticlesandDepose (int lev, amrex::Real cur_time); void PushParticlesandDepose ( amrex::Real cur_time); @@ -486,6 +490,8 @@ private: int do_pml = 1; int pml_ncell = 10; int pml_delta = 10; + int pml_has_particles = 0; + int do_pml_j_damping = 0; amrex::Vector<std::unique_ptr<PML> > pml; amrex::Real moving_window_x = std::numeric_limits<amrex::Real>::max(); @@ -496,7 +502,7 @@ private: int warpx_do_continuous_injection = 0; int num_injected_species = -1; amrex::Vector<int> injected_plasma_species; - + int do_electrostatic = 0; int n_buffer = 4; amrex::Real const_dt = 0.5e-11; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 877882037..77fd83f7a 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -380,6 +380,8 @@ WarpX::ReadParameters () pp.query("do_pml", do_pml); pp.query("pml_ncell", pml_ncell); pp.query("pml_delta", pml_delta); + pp.query("pml_has_particles", pml_has_particles); + pp.query("do_pml_j_damping", do_pml_j_damping); pp.query("dump_openpmd", dump_openpmd); pp.query("dump_plotfiles", dump_plotfiles); |