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