diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 0cea2709312..281aa5e75ba 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -57,6 +57,199 @@ using namespace amrex; +namespace +{ + + /** This function shifts a MultiFab in a given direction + * + * \param[in,out] mf the MultiFab to be shifted + * \param[in] geom the Geometry object associated to the level of the MultiFab mf + * \param[in] num_shift magnitude of the shift (cell number) + * \param[in] dir direction of the shift + * \param[in] safe_guard_cells flag to enable "safe mode" data exchanges with more guard cells + * \param[in] do_single_precision_comms flag to enable single precision communications + * \param[in,out] cost the pointer to the data structure holding costs for timer-based load-balance + * \param[in] external_field the external field (used to initialize EM fields) + * \param[in] useparser flag to enable the use of a field parser to initialize EM fields + * \param[in] field_parser the field parser + * \param[in] PMLRZ_flag flag to enable a special treatment for PML in RZ simulations + */ + void shiftMF ( + amrex::MultiFab& mf, const amrex::Geometry& geom, + int num_shift, int dir, + bool safe_guard_cells, bool do_single_precision_comms, + amrex::LayoutData* cost, + amrex::Real external_field=0.0, bool useparser = false, + amrex::ParserExecutor<3> const& field_parser={}, + const bool PMLRZ_flag = false) + { + using namespace amrex::literals; + WARPX_PROFILE("warpx::shiftMF()"); + const amrex::BoxArray& ba = mf.boxArray(); + const amrex::DistributionMapping& dm = mf.DistributionMap(); + const int nc = mf.nComp(); + const amrex::IntVect& ng = mf.nGrowVect(); + + AMREX_ALWAYS_ASSERT(ng[dir] >= std::abs(num_shift)); + + amrex::MultiFab tmpmf(ba, dm, nc, ng); + amrex::MultiFab::Copy(tmpmf, mf, 0, 0, nc, ng); + + if ( safe_guard_cells ) { + // Fill guard cells. + ablastr::utils::communication::FillBoundary(tmpmf, do_single_precision_comms, geom.periodicity()); + } else { + amrex::IntVect ng_mw = amrex::IntVect::TheUnitVector(); + // Enough guard cells in the MW direction + ng_mw[dir] = std::abs(num_shift); + // Make sure we don't exceed number of guard cells allocated + ng_mw = ng_mw.min(ng); + // Fill guard cells. + ablastr::utils::communication::FillBoundary(tmpmf, ng_mw, do_single_precision_comms, geom.periodicity()); + } + + // Make a box that covers the region that the window moved into + const amrex::IndexType& typ = ba.ixType(); + const amrex::Box& domainBox = geom.Domain(); + amrex::Box adjBox; + if (num_shift > 0) { + adjBox = adjCellHi(domainBox, dir, ng[dir]); + } else { + adjBox = adjCellLo(domainBox, dir, ng[dir]); + } + adjBox = amrex::convert(adjBox, typ); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if (idim == dir and typ.nodeCentered(dir)) { + if (num_shift > 0) { + adjBox.growLo(idim, -1); + } else { + adjBox.growHi(idim, -1); + } + } else if (idim != dir) { + adjBox.growLo(idim, ng[idim]); + adjBox.growHi(idim, ng[idim]); + } + } + + amrex::IntVect shiftiv(0); + shiftiv[dir] = num_shift; + const amrex::Dim3 shift = shiftiv.dim3(); + + const amrex::RealBox& real_box = geom.ProbDomain(); + const auto dx = geom.CellSizeArray(); + +#ifdef AMREX_USE_OMP + #pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (amrex::MFIter mfi(tmpmf, TilingIfNotGPU()); mfi.isValid(); ++mfi ) + { + if (cost) + { + amrex::Gpu::synchronize(); + } + auto wt = static_cast(amrex::second()); + + auto const& dstfab = mf.array(mfi); + auto const& srcfab = tmpmf.array(mfi); + + const amrex::Box& outbox = mfi.growntilebox() & adjBox; + + if (outbox.ok()) { + if (!useparser) { + AMREX_PARALLEL_FOR_4D ( outbox, nc, i, j, k, n, + { + srcfab(i,j,k,n) = external_field; + }) + } else { + // index type of the src mf + auto const& mf_IndexType = (tmpmf).ixType(); + amrex::IntVect mf_type(AMREX_D_DECL(0,0,0)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + mf_type[idim] = mf_IndexType.nodeCentered(idim); + } + + amrex::ParallelFor (outbox, nc, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + // Compute x,y,z co-ordinates based on index type of mf +#if defined(WARPX_DIM_1D_Z) + const amrex::Real x = 0.0_rt; + const amrex::Real y = 0.0_rt; + const amrex::Real fac_z = (1.0_rt - mf_type[0]) * dx[0]*0.5_rt; + const amrex::Real z = i*dx[0] + real_box.lo(0) + fac_z; +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) + const amrex::Real fac_x = (1.0_rt - mf_type[0]) * dx[0]*0.5_rt; + const amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; + const amrex::Real y = 0.0; + const amrex::Real fac_z = (1.0_rt - mf_type[1]) * dx[1]*0.5_rt; + const amrex::Real z = j*dx[1] + real_box.lo(1) + fac_z; +#else + const amrex::Real fac_x = (1.0_rt - mf_type[0]) * dx[0]*0.5_rt; + const amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; + const amrex::Real fac_y = (1.0_rt - mf_type[1]) * dx[1]*0.5_rt; + const amrex::Real y = j*dx[1] + real_box.lo(1) + fac_y; + const amrex::Real fac_z = (1.0_rt - mf_type[2]) * dx[2]*0.5_rt; + const amrex::Real z = k*dx[2] + real_box.lo(2) + fac_z; +#endif + srcfab(i,j,k,n) = field_parser(x,y,z); + }); + } + + } + + amrex::Box dstBox = mf[mfi].box(); + if (num_shift > 0) { + dstBox.growHi(dir, -num_shift); + } else { + dstBox.growLo(dir, num_shift); + } + AMREX_PARALLEL_FOR_4D ( dstBox, nc, i, j, k, n, + { + dstfab(i,j,k,n) = srcfab(i+shift.x,j+shift.y,k+shift.z,n); + }) + + if (cost) + { + amrex::Gpu::synchronize(); + wt = static_cast(amrex::second()) - wt; + amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt); + } + } + +#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) + if (PMLRZ_flag) { + // This does the exchange of data in the corner guard cells, the cells that are in the + // guard region both radially and longitudinally. These are the PML cells in the overlapping + // longitudinal region. FillBoundary normally does not update these cells. + // This update is needed so that the cells at the end of the FABs are updated appropriately + // with the data shifted from the neighboring FAB. Without this update, the RZ PML becomes + // unstable with the moving grid. + // This code creates a temporary MultiFab using a BoxList where the radial size of all of + // its boxes is increased so that the radial guard cells are included in the boxes valid domain. + // The temporary MultiFab is setup to refer to the data of the original Multifab (this can + // be done since the shape of the data is all the same, just the indexing is different). + amrex::BoxList bl; + const auto ba_size = static_cast(ba.size()); + for (int i = 0; i < ba_size; ++i) { + bl.push_back(amrex::grow(ba[i], 0, mf.nGrowVect()[0])); + } + const amrex::BoxArray rba(std::move(bl)); + amrex::MultiFab rmf(rba, dm, mf.nComp(), IntVect(0,mf.nGrowVect()[1]), MFInfo().SetAlloc(false)); + + for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) { + rmf.setFab(mfi, FArrayBox(mf[mfi], amrex::make_alias, 0, mf.nComp())); + } + rmf.FillBoundary(false); + } +#else + amrex::ignore_unused(PMLRZ_flag); +#endif + + } + +} + void WarpX::UpdateInjectionPosition (const amrex::Real a_dt) { @@ -208,9 +401,6 @@ WarpX::MoveWindow (const int step, bool move_j) int num_shift = num_shift_base; int num_shift_crse = num_shift; - constexpr auto do_update_cost = true; - constexpr auto dont_update_cost = false; //We can't update cost for PML - // Shift the mesh fields for (int lev = 0; lev <= finest_level; ++lev) { @@ -219,6 +409,11 @@ WarpX::MoveWindow (const int step, bool move_j) num_shift *= refRatio(lev-1)[dir]; } + auto* cost_lev = + (WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) ? getCosts(lev) : nullptr; + + amrex::LayoutData* no_cost = nullptr ; //We can't update cost for PML + // Shift each component of vector fields (E, B, j) for (int dim = 0; dim < 3; ++dim) { // Fine grid @@ -240,59 +435,60 @@ WarpX::MoveWindow (const int step, bool move_j) if (dim == 1) { Efield_parser = m_p_ext_field_params->Eyfield_parser->compile<3>(); } if (dim == 2) { Efield_parser = m_p_ext_field_params->Ezfield_parser->compile<3>(); } } - shiftMF(*m_fields.get(FieldType::Bfield_fp, Direction{dim}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*m_fields.get(FieldType::Bfield_fp, Direction{dim}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params->B_external_grid[dim], use_Bparser, Bfield_parser); - shiftMF(*m_fields.get(FieldType::Efield_fp, Direction{dim}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*m_fields.get(FieldType::Efield_fp, Direction{dim}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params->E_external_grid[dim], use_Eparser, Efield_parser); if (fft_do_time_averaging) { ablastr::fields::MultiLevelVectorField Efield_avg_fp = m_fields.get_mr_levels_alldirs(FieldType::Efield_avg_fp, finest_level); ablastr::fields::MultiLevelVectorField Bfield_avg_fp = m_fields.get_mr_levels_alldirs(FieldType::Bfield_avg_fp, finest_level); - shiftMF(*Bfield_avg_fp[lev][dim], geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*Bfield_avg_fp[lev][dim], geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params->B_external_grid[dim], use_Bparser, Bfield_parser); - shiftMF(*Efield_avg_fp[lev][dim], geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*Efield_avg_fp[lev][dim], geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params-> E_external_grid[dim], use_Eparser, Efield_parser); } if (move_j) { - shiftMF(*m_fields.get(FieldType::current_fp, Direction{dim}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::current_fp, Direction{dim}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); } if (pml[lev] && pml[lev]->ok()) { amrex::MultiFab* pml_B = m_fields.get(FieldType::pml_B_fp, Direction{dim}, lev); amrex::MultiFab* pml_E = m_fields.get(FieldType::pml_E_fp, Direction{dim}, lev); - shiftMF(*pml_B, geom[lev], num_shift, dir, lev, dont_update_cost, m_safe_guard_cells); - shiftMF(*pml_E, geom[lev], num_shift, dir, lev, dont_update_cost, m_safe_guard_cells); + ::shiftMF(*pml_B, geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); + ::shiftMF(*pml_E, geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); } #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) + const bool PMLRZ_flag = getPMLRZ(); if (pml_rz[lev] && dim < 2) { amrex::MultiFab* pml_rz_B = m_fields.get(FieldType::pml_B_fp, Direction{dim}, lev); amrex::MultiFab* pml_rz_E = m_fields.get(FieldType::pml_E_fp, Direction{dim}, lev); - shiftMF(*pml_rz_B, geom[lev], num_shift, dir, lev, dont_update_cost, m_safe_guard_cells); - shiftMF(*pml_rz_E, geom[lev], num_shift, dir, lev, dont_update_cost, m_safe_guard_cells); + ::shiftMF(*pml_rz_B, geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, no_cost, 0.0_rt, false, amrex::ParserExecutor<3>{}, PMLRZ_flag); + ::shiftMF(*pml_rz_E, geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, no_cost, 0.0_rt, false, amrex::ParserExecutor<3>{}, PMLRZ_flag); } #endif if (lev > 0) { // coarse grid - shiftMF(*m_fields.get(FieldType::Bfield_cp, Direction{dim}, lev), geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*m_fields.get(FieldType::Bfield_cp, Direction{dim}, lev), geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params->B_external_grid[dim], use_Bparser, Bfield_parser); - shiftMF(*m_fields.get(FieldType::Efield_cp, Direction{dim}, lev), geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*m_fields.get(FieldType::Efield_cp, Direction{dim}, lev), geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params->E_external_grid[dim], use_Eparser, Efield_parser); - shiftMF(*m_fields.get(FieldType::Bfield_aux, Direction{dim}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells); - shiftMF(*m_fields.get(FieldType::Efield_aux, Direction{dim}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::Bfield_aux, Direction{dim}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); + ::shiftMF(*m_fields.get(FieldType::Efield_aux, Direction{dim}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); if (fft_do_time_averaging) { ablastr::fields::MultiLevelVectorField Efield_avg_cp = m_fields.get_mr_levels_alldirs(FieldType::Efield_avg_cp, finest_level, skip_lev0_coarse_patch); ablastr::fields::MultiLevelVectorField Bfield_avg_cp = m_fields.get_mr_levels_alldirs(FieldType::Bfield_avg_cp, finest_level, skip_lev0_coarse_patch); - shiftMF(*Bfield_avg_cp[lev][dim], geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*Bfield_avg_cp[lev][dim], geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params->B_external_grid[dim], use_Bparser, Bfield_parser); - shiftMF(*Efield_avg_cp[lev][dim], geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells, + ::shiftMF(*Efield_avg_cp[lev][dim], geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev, m_p_ext_field_params->E_external_grid[dim], use_Eparser, Efield_parser); } if (move_j) { - shiftMF(*m_fields.get(FieldType::current_cp, Direction{dim}, lev), geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::current_cp, Direction{dim}, lev), geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); } if (do_pml && pml[lev]->ok()) { amrex::MultiFab* pml_B_cp = m_fields.get(FieldType::pml_B_cp, Direction{dim}, lev); amrex::MultiFab* pml_E_cp = m_fields.get(FieldType::pml_E_cp, Direction{dim}, lev); - shiftMF(*pml_B_cp, geom[lev-1], num_shift_crse, dir, lev, dont_update_cost, m_safe_guard_cells); - shiftMF(*pml_E_cp, geom[lev-1], num_shift_crse, dir, lev, dont_update_cost, m_safe_guard_cells); + ::shiftMF(*pml_B_cp, geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); + ::shiftMF(*pml_E_cp, geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); } } } @@ -302,11 +498,11 @@ WarpX::MoveWindow (const int step, bool move_j) if (m_fields.has(FieldType::F_fp, lev)) { // Fine grid - shiftMF(*m_fields.get(FieldType::F_fp, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::F_fp, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); if (lev > 0) { // Coarse grid - shiftMF(*m_fields.get(FieldType::F_cp, lev), geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::F_cp, lev), geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); } } @@ -317,7 +513,7 @@ WarpX::MoveWindow (const int step, bool move_j) if (do_pml && pml[lev]->ok()) { amrex::MultiFab* pml_F = m_fields.get(FieldType::pml_F_fp, lev); - shiftMF(*pml_F, geom[lev], num_shift, dir, lev, dont_update_cost, m_safe_guard_cells); + ::shiftMF(*pml_F, geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); } if (lev > 0) { @@ -325,7 +521,7 @@ WarpX::MoveWindow (const int step, bool move_j) if (do_pml && pml[lev]->ok()) { amrex::MultiFab* pml_F = m_fields.get(FieldType::pml_F_cp, lev); - shiftMF(*pml_F, geom[lev-1], num_shift_crse, dir, lev, dont_update_cost, m_safe_guard_cells); + ::shiftMF(*pml_F, geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); } } } @@ -335,11 +531,11 @@ WarpX::MoveWindow (const int step, bool move_j) if (m_fields.has(FieldType::G_fp, lev)) { // Fine grid - shiftMF(*m_fields.get(FieldType::G_fp, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::G_fp, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); if (lev > 0) { // Coarse grid - shiftMF(*m_fields.get(FieldType::G_cp, lev), geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::G_cp, lev), geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); } } @@ -350,7 +546,7 @@ WarpX::MoveWindow (const int step, bool move_j) if (do_pml && pml[lev]->ok()) { amrex::MultiFab* pml_G = m_fields.get(FieldType::pml_G_fp, lev); - shiftMF(*pml_G, geom[lev], num_shift, dir, lev, dont_update_cost, m_safe_guard_cells); + ::shiftMF(*pml_G, geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); } if (lev > 0) { @@ -358,7 +554,7 @@ WarpX::MoveWindow (const int step, bool move_j) if (do_pml && pml[lev]->ok()) { amrex::MultiFab* pml_G = m_fields.get(FieldType::pml_G_cp, lev); - shiftMF(*pml_G, geom[lev-1], num_shift_crse, dir, lev, dont_update_cost, m_safe_guard_cells); + ::shiftMF(*pml_G, geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, no_cost); } } } @@ -367,10 +563,10 @@ WarpX::MoveWindow (const int step, bool move_j) if (move_j) { if (m_fields.has(FieldType::rho_fp, lev)) { // Fine grid - shiftMF(*m_fields.get(FieldType::rho_fp,lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::rho_fp,lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); if (lev > 0){ // Coarse grid - shiftMF(*m_fields.get(FieldType::rho_cp,lev), geom[lev-1], num_shift_crse, dir, lev, do_update_cost, m_safe_guard_cells); + ::shiftMF(*m_fields.get(FieldType::rho_cp,lev), geom[lev-1], num_shift_crse, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev); } } } @@ -380,10 +576,10 @@ WarpX::MoveWindow (const int step, bool move_j) const int n_fluid_species = myfl->nSpecies(); for (int i=0; iGetFluidContainer(i); - shiftMF( *m_fields.get(fl.name_mf_N, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells ); - shiftMF( *m_fields.get(fl.name_mf_NU, Direction{0}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells ); - shiftMF( *m_fields.get(fl.name_mf_NU, Direction{1}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells ); - shiftMF( *m_fields.get(fl.name_mf_NU, Direction{2}, lev), geom[lev], num_shift, dir, lev, do_update_cost, m_safe_guard_cells ); + ::shiftMF( *m_fields.get(fl.name_mf_N, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev ); + ::shiftMF( *m_fields.get(fl.name_mf_NU, Direction{0}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev ); + ::shiftMF( *m_fields.get(fl.name_mf_NU, Direction{1}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev ); + ::shiftMF( *m_fields.get(fl.name_mf_NU, Direction{2}, lev), geom[lev], num_shift, dir, m_safe_guard_cells, do_single_precision_comms, cost_lev ); } } } @@ -477,179 +673,6 @@ WarpX::MoveWindow (const int step, bool move_j) return num_shift_base; } -void -WarpX::shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, - int num_shift, int dir, const int lev, bool update_cost_flag, - const bool safe_guard_cells, - amrex::Real external_field, bool useparser, - amrex::ParserExecutor<3> const& field_parser) -{ - using namespace amrex::literals; - WARPX_PROFILE("WarpX::shiftMF()"); - const amrex::BoxArray& ba = mf.boxArray(); - const amrex::DistributionMapping& dm = mf.DistributionMap(); - const int nc = mf.nComp(); - const amrex::IntVect& ng = mf.nGrowVect(); - - AMREX_ALWAYS_ASSERT(ng[dir] >= num_shift); - - amrex::MultiFab tmpmf(ba, dm, nc, ng); - amrex::MultiFab::Copy(tmpmf, mf, 0, 0, nc, ng); - - if ( safe_guard_cells ) { - // Fill guard cells. - ablastr::utils::communication::FillBoundary(tmpmf, WarpX::do_single_precision_comms, geom.periodicity()); - } else { - amrex::IntVect ng_mw = amrex::IntVect::TheUnitVector(); - // Enough guard cells in the MW direction - ng_mw[dir] = num_shift; - // Make sure we don't exceed number of guard cells allocated - ng_mw = ng_mw.min(ng); - // Fill guard cells. - ablastr::utils::communication::FillBoundary(tmpmf, ng_mw, WarpX::do_single_precision_comms, geom.periodicity()); - } - - // Make a box that covers the region that the window moved into - const amrex::IndexType& typ = ba.ixType(); - const amrex::Box& domainBox = geom.Domain(); - amrex::Box adjBox; - if (num_shift > 0) { - adjBox = adjCellHi(domainBox, dir, ng[dir]); - } else { - adjBox = adjCellLo(domainBox, dir, ng[dir]); - } - adjBox = amrex::convert(adjBox, typ); - - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if (idim == dir and typ.nodeCentered(dir)) { - if (num_shift > 0) { - adjBox.growLo(idim, -1); - } else { - adjBox.growHi(idim, -1); - } - } else if (idim != dir) { - adjBox.growLo(idim, ng[idim]); - adjBox.growHi(idim, ng[idim]); - } - } - - amrex::IntVect shiftiv(0); - shiftiv[dir] = num_shift; - const amrex::Dim3 shift = shiftiv.dim3(); - - const amrex::RealBox& real_box = geom.ProbDomain(); - const auto dx = geom.CellSizeArray(); - - amrex::LayoutData* cost = WarpX::getCosts(lev); -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - - for (amrex::MFIter mfi(tmpmf, TilingIfNotGPU()); mfi.isValid(); ++mfi ) - { - if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) - { - amrex::Gpu::synchronize(); - } - auto wt = static_cast(amrex::second()); - - auto const& dstfab = mf.array(mfi); - auto const& srcfab = tmpmf.array(mfi); - - const amrex::Box& outbox = mfi.growntilebox() & adjBox; - - if (outbox.ok()) { - if (!useparser) { - AMREX_PARALLEL_FOR_4D ( outbox, nc, i, j, k, n, - { - srcfab(i,j,k,n) = external_field; - }) - } else { - // index type of the src mf - auto const& mf_IndexType = (tmpmf).ixType(); - amrex::IntVect mf_type(AMREX_D_DECL(0,0,0)); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - mf_type[idim] = mf_IndexType.nodeCentered(idim); - } - - amrex::ParallelFor (outbox, nc, - [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - // Compute x,y,z co-ordinates based on index type of mf -#if defined(WARPX_DIM_1D_Z) - const amrex::Real x = 0.0_rt; - const amrex::Real y = 0.0_rt; - const amrex::Real fac_z = (1.0_rt - mf_type[0]) * dx[0]*0.5_rt; - const amrex::Real z = i*dx[0] + real_box.lo(0) + fac_z; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const amrex::Real fac_x = (1.0_rt - mf_type[0]) * dx[0]*0.5_rt; - const amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; - const amrex::Real y = 0.0; - const amrex::Real fac_z = (1.0_rt - mf_type[1]) * dx[1]*0.5_rt; - const amrex::Real z = j*dx[1] + real_box.lo(1) + fac_z; -#else - const amrex::Real fac_x = (1.0_rt - mf_type[0]) * dx[0]*0.5_rt; - const amrex::Real x = i*dx[0] + real_box.lo(0) + fac_x; - const amrex::Real fac_y = (1.0_rt - mf_type[1]) * dx[1]*0.5_rt; - const amrex::Real y = j*dx[1] + real_box.lo(1) + fac_y; - const amrex::Real fac_z = (1.0_rt - mf_type[2]) * dx[2]*0.5_rt; - const amrex::Real z = k*dx[2] + real_box.lo(2) + fac_z; -#endif - srcfab(i,j,k,n) = field_parser(x,y,z); - }); - } - - } - - amrex::Box dstBox = mf[mfi].box(); - if (num_shift > 0) { - dstBox.growHi(dir, -num_shift); - } else { - dstBox.growLo(dir, num_shift); - } - AMREX_PARALLEL_FOR_4D ( dstBox, nc, i, j, k, n, - { - dstfab(i,j,k,n) = srcfab(i+shift.x,j+shift.y,k+shift.z,n); - }) - - if (cost && update_cost_flag && - WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers) - { - amrex::Gpu::synchronize(); - wt = static_cast(amrex::second()) - wt; - amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt); - } - } - -#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_FFT) - if (WarpX::GetInstance().getPMLRZ()) { - // This does the exchange of data in the corner guard cells, the cells that are in the - // guard region both radially and longitudinally. These are the PML cells in the overlapping - // longitudinal region. FillBoundary normally does not update these cells. - // This update is needed so that the cells at the end of the FABs are updated appropriately - // with the data shifted from the neighboring FAB. Without this update, the RZ PML becomes - // unstable with the moving grid. - // This code creates a temporary MultiFab using a BoxList where the radial size of all of - // its boxes is increased so that the radial guard cells are included in the boxes valid domain. - // The temporary MultiFab is setup to refer to the data of the original Multifab (this can - // be done since the shape of the data is all the same, just the indexing is different). - amrex::BoxList bl; - const auto ba_size = static_cast(ba.size()); - for (int i = 0; i < ba_size; ++i) { - bl.push_back(amrex::grow(ba[i], 0, mf.nGrowVect()[0])); - } - const amrex::BoxArray rba(std::move(bl)); - amrex::MultiFab rmf(rba, dm, mf.nComp(), IntVect(0,mf.nGrowVect()[1]), MFInfo().SetAlloc(false)); - - for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) { - rmf.setFab(mfi, FArrayBox(mf[mfi], amrex::make_alias, 0, mf.nComp())); - } - rmf.FillBoundary(false); - } -#endif - -} - void WarpX::ShiftGalileanBoundary () { diff --git a/Source/WarpX.H b/Source/WarpX.H index b12cb1ab7f0..7d164a9e685 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -166,12 +166,6 @@ public: amrex::Vector,3 > >& GetEBUpdateEFlag() { return m_eb_update_E; } amrex::Vector< std::unique_ptr > const & GetEBReduceParticleShapeFlag() const { return m_eb_reduce_particle_shape; } - static void shiftMF (amrex::MultiFab& mf, const amrex::Geometry& geom, - int num_shift, int dir, int lev, bool update_cost_flag, - bool safe_guard_cells, - amrex::Real external_field=0.0, bool useparser = false, - amrex::ParserExecutor<3> const& field_parser={}); - /** * \brief * If an authors' string is specified in the inputfile, this method returns that string.