diff options
Diffstat (limited to 'Source/Parallelization')
-rw-r--r-- | Source/Parallelization/WarpXComm.cpp | 151 | ||||
-rw-r--r-- | Source/Parallelization/WarpXComm_K.H | 485 | ||||
-rw-r--r-- | Source/Parallelization/WarpXRegrid.cpp | 2 |
3 files changed, 637 insertions, 1 deletions
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 92f0b4f09..52df3dc25 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -54,6 +54,157 @@ WarpX::UpdateAuxilaryData () { BL_PROFILE("UpdateAuxilaryData()"); + if (Bfield_aux[0][0]->ixType() == Bfield_fp[0][0]->ixType()) { + UpdateAuxilaryDataSameType(); + } else { + UpdateAuxilaryDataStagToNodal(); + } +} + +void +WarpX::UpdateAuxilaryDataStagToNodal () +{ + // For level 0, we only need to do the average. +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Bfield_aux[0][0]); mfi.isValid(); ++mfi) + { + Array4<Real> const& bx_aux = Bfield_aux[0][0]->array(mfi); + Array4<Real> const& by_aux = Bfield_aux[0][1]->array(mfi); + Array4<Real> const& bz_aux = Bfield_aux[0][2]->array(mfi); + Array4<Real const> const& bx_fp = Bfield_fp[0][0]->const_array(mfi); + Array4<Real const> const& by_fp = Bfield_fp[0][1]->const_array(mfi); + Array4<Real const> const& bz_fp = Bfield_fp[0][2]->const_array(mfi); + + Array4<Real> const& ex_aux = Efield_aux[0][0]->array(mfi); + Array4<Real> const& ey_aux = Efield_aux[0][1]->array(mfi); + Array4<Real> const& ez_aux = Efield_aux[0][2]->array(mfi); + Array4<Real const> const& ex_fp = Efield_fp[0][0]->const_array(mfi); + Array4<Real const> const& ey_fp = Efield_fp[0][1]->const_array(mfi); + Array4<Real const> const& ez_fp = Efield_fp[0][2]->const_array(mfi); + + const Box& bx = mfi.fabbox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp); + warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp); + warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp); + warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp); + warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp); + warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp); + }); + } + + for (int lev = 1; lev <= finest_level; ++lev) + { + BoxArray const& nba = Bfield_aux[lev][0]->boxArray(); + BoxArray const& cnba = amrex::coarsen(nba,2); + DistributionMapping const& dm = Bfield_aux[lev][0]->DistributionMap(); + auto const& cperiod = Geom(lev-1).periodicity(); + + // Bfield + { + Array<std::unique_ptr<MultiFab>,3> Btmp; + if (Bfield_cax[lev][0]) { + for (int i = 0; i < 3; ++i) { + Btmp[i].reset(new MultiFab(*Bfield_cax[lev][i], amrex::make_alias, 0, 1)); + } + } else { + IntVect ngtmp = Bfield_aux[lev-1][0]->nGrowVect(); + for (int i = 0; i < 3; ++i) { + Btmp[i].reset(new MultiFab(cnba, dm, 1, ngtmp)); + } + } + // ParallelCopy from coarse level + for (int i = 0; i < 3; ++i) { + IntVect ng = Btmp[i]->nGrowVect(); + Btmp[i]->ParallelCopy(*Bfield_aux[lev-1][i], 0, 0, 1, ng, ng, cperiod); + } + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Bfield_aux[lev][0]); mfi.isValid(); ++mfi) + { + Array4<Real> const& bx_aux = Bfield_aux[lev][0]->array(mfi); + Array4<Real> const& by_aux = Bfield_aux[lev][1]->array(mfi); + Array4<Real> const& bz_aux = Bfield_aux[lev][2]->array(mfi); + Array4<Real const> const& bx_fp = Bfield_fp[lev][0]->const_array(mfi); + Array4<Real const> const& by_fp = Bfield_fp[lev][1]->const_array(mfi); + Array4<Real const> const& bz_fp = Bfield_fp[lev][2]->const_array(mfi); + Array4<Real const> const& bx_cp = Bfield_cp[lev][0]->const_array(mfi); + Array4<Real const> const& by_cp = Bfield_cp[lev][1]->const_array(mfi); + Array4<Real const> const& bz_cp = Bfield_cp[lev][2]->const_array(mfi); + Array4<Real const> const& bx_c = Btmp[0]->const_array(mfi); + Array4<Real const> const& by_c = Btmp[1]->const_array(mfi); + Array4<Real const> const& bz_c = Btmp[2]->const_array(mfi); + + const Box& bx = mfi.fabbox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp, bx_cp, bx_c); + warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp, by_cp, by_c); + warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp, bz_cp, bz_c); + }); + } + } + + // Efield + { + Array<std::unique_ptr<MultiFab>,3> Etmp; + if (Efield_cax[lev][0]) { + for (int i = 0; i < 3; ++i) { + Etmp[i].reset(new MultiFab(*Efield_cax[lev][i], amrex::make_alias, 0, 1)); + } + } else { + IntVect ngtmp = Efield_aux[lev-1][0]->nGrowVect(); + for (int i = 0; i < 3; ++i) { + Etmp[i].reset(new MultiFab(cnba, dm, 1, ngtmp)); + } + } + // ParallelCopy from coarse level + for (int i = 0; i < 3; ++i) { + IntVect ng = Etmp[i]->nGrowVect(); + Etmp[i]->ParallelCopy(*Efield_aux[lev-1][i], 0, 0, 1, ng, ng, cperiod); + } + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi) + { + Array4<Real> const& ex_aux = Efield_aux[lev][0]->array(mfi); + Array4<Real> const& ey_aux = Efield_aux[lev][1]->array(mfi); + Array4<Real> const& ez_aux = Efield_aux[lev][2]->array(mfi); + Array4<Real const> const& ex_fp = Efield_fp[lev][0]->const_array(mfi); + Array4<Real const> const& ey_fp = Efield_fp[lev][1]->const_array(mfi); + Array4<Real const> const& ez_fp = Efield_fp[lev][2]->const_array(mfi); + Array4<Real const> const& ex_cp = Efield_cp[lev][0]->const_array(mfi); + Array4<Real const> const& ey_cp = Efield_cp[lev][1]->const_array(mfi); + Array4<Real const> const& ez_cp = Efield_cp[lev][2]->const_array(mfi); + Array4<Real const> const& ex_c = Etmp[0]->const_array(mfi); + Array4<Real const> const& ey_c = Etmp[1]->const_array(mfi); + Array4<Real const> const& ez_c = Etmp[2]->const_array(mfi); + + const Box& bx = mfi.fabbox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp, ex_cp, ex_c); + warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp, ey_cp, ey_c); + warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp, ez_cp, ez_c); + }); + } + } + } +} + +void +WarpX::UpdateAuxilaryDataSameType () +{ for (int lev = 1; lev <= finest_level; ++lev) { const auto& crse_period = Geom(lev-1).periodicity(); diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index 093323ec3..5da867c9f 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -158,4 +158,489 @@ void warpx_interp_efield_z (int j, int k, int l, #endif } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bxa, + amrex::Array4<amrex::Real const> const& Bxf, + amrex::Array4<amrex::Real const> const& Bxc, + amrex::Array4<amrex::Real const> const& Bxg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * Bxg(jg ,kg ,0) + + owx * wy * Bxg(jg ,kg+1,0) + + wx * owy * Bxg(jg+1,kg ,0) + + wx * wy * Bxg(jg+1,kg+1,0); + + // interp from coarse staggered to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + Real bc = owx * owy * Bxc(jg ,kg ,0) + + owx * wy * Bxc(jg ,kg-1,0) + + wx * owy * Bxc(jg+1,kg ,0) + + wx * wy * Bxc(jg+1,kg-1,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5*(Bxf(j,k-1,0) + Bxf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * owz * Bxg(jg ,kg ,lg ) + + wx * owy * owz * Bxg(jg+1,kg ,lg ) + + owx * wy * owz * Bxg(jg ,kg+1,lg ) + + wx * wy * owz * Bxg(jg+1,kg+1,lg ) + + owx * owy * wz * Bxg(jg ,kg ,lg+1) + + wx * owy * wz * Bxg(jg+1,kg ,lg+1) + + owx * wy * wz * Bxg(jg ,kg+1,lg+1) + + wx * wy * wz * Bxg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + wz = 0.5-wz; owz = 1.0-wz; + Real bc = owx * owy * owz * Bxc(jg ,kg ,lg ) + + wx * owy * owz * Bxc(jg+1,kg ,lg ) + + owx * wy * owz * Bxc(jg ,kg-1,lg ) + + wx * wy * owz * Bxc(jg+1,kg-1,lg ) + + owx * owy * wz * Bxc(jg ,kg ,lg-1) + + wx * owy * wz * Bxc(jg+1,kg ,lg-1) + + owx * wy * wz * Bxc(jg ,kg-1,lg-1) + + wx * wy * wz * Bxc(jg+1,kg-1,lg-1); + + // interp from fine stagged to fine nodal + Real bf = 0.25*(Bxf(j,k-1,l-1) + Bxf(j,k,l-1) + Bxf(j,k-1,l) + Bxf(j,k,l)); +#endif + + Bxa(j,k,l) = bg + (bf-bc); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bya, + amrex::Array4<amrex::Real const> const& Byf, + amrex::Array4<amrex::Real const> const& Byc, + amrex::Array4<amrex::Real const> const& Byg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * Byg(jg ,kg ,0) + + owx * wy * Byg(jg ,kg+1,0) + + wx * owy * Byg(jg+1,kg ,0) + + wx * wy * Byg(jg+1,kg+1,0); + + // interp from coarse stagged (cell-centered for By) to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + wy = 0.5-wy; owy = 1.0-wy; + Real bc = owx * owy * Byc(jg ,kg ,0) + + owx * wy * Byc(jg ,kg-1,0) + + wx * owy * Byc(jg-1,kg ,0) + + wx * wy * Byc(jg-1,kg-1,0); + + // interp form fine stagger (cell-centered for By) to fine nodal + Real bf = 0.25*(Byf(j,k,0) + Byf(j-1,k,0) + Byf(j,k-1,0) + Byf(j-1,k-1,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * owz * Byg(jg ,kg ,lg ) + + wx * owy * owz * Byg(jg+1,kg ,lg ) + + owx * wy * owz * Byg(jg ,kg+1,lg ) + + wx * wy * owz * Byg(jg+1,kg+1,lg ) + + owx * owy * wz * Byg(jg ,kg ,lg+1) + + wx * owy * wz * Byg(jg+1,kg ,lg+1) + + owx * wy * wz * Byg(jg ,kg+1,lg+1) + + wx * wy * wz * Byg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + wz = 0.5-wz; owz = 1.0-wz; + Real bc = owx * owy * owz * Byc(jg ,kg ,lg ) + + wx * owy * owz * Byc(jg-1,kg ,lg ) + + owx * wy * owz * Byc(jg ,kg+1,lg ) + + wx * wy * owz * Byc(jg-1,kg+1,lg ) + + owx * owy * wz * Byc(jg ,kg ,lg-1) + + wx * owy * wz * Byc(jg-1,kg ,lg-1) + + owx * wy * wz * Byc(jg ,kg+1,lg-1) + + wx * wy * wz * Byc(jg-1,kg+1,lg-1); + + // interp from fine stagged to fine nodal + Real bf = 0.25*(Byf(j-1,k,l-1) + Byf(j,k,l-1) + Byf(j-1,k,l) + Byf(j,k,l)); + +#endif + + Bya(j,k,l) = bg + (bf-bc); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bza, + amrex::Array4<amrex::Real const> const& Bzf, + amrex::Array4<amrex::Real const> const& Bzc, + amrex::Array4<amrex::Real const> const& Bzg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * Bzg(jg ,kg ,0) + + owx * wy * Bzg(jg ,kg+1,0) + + wx * owy * Bzg(jg+1,kg ,0) + + wx * wy * Bzg(jg+1,kg+1,0); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + Real bc = owx * owy * Bzc(jg ,kg ,0) + + owx * wy * Bzc(jg ,kg+1,0) + + wx * owy * Bzc(jg-1,kg ,0) + + wx * wy * Bzc(jg-1,kg+1,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5*(Bzf(j-1,k,0) + Bzf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * owz * Bzg(jg ,kg ,lg ) + + wx * owy * owz * Bzg(jg+1,kg ,lg ) + + owx * wy * owz * Bzg(jg ,kg+1,lg ) + + wx * wy * owz * Bzg(jg+1,kg+1,lg ) + + owx * owy * wz * Bzg(jg ,kg ,lg+1) + + wx * owy * wz * Bzg(jg+1,kg ,lg+1) + + owx * wy * wz * Bzg(jg ,kg+1,lg+1) + + wx * wy * wz * Bzg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + wy = 0.5-wy; owy = 1.0-wy; + Real bc = owx * owy * owz * Bzc(jg ,kg ,lg ) + + wx * owy * owz * Bzc(jg-1,kg ,lg ) + + owx * wy * owz * Bzc(jg ,kg-1,lg ) + + wx * wy * owz * Bzc(jg-1,kg-1,lg ) + + owx * owy * wz * Bzc(jg ,kg ,lg+1) + + wx * owy * wz * Bzc(jg-1,kg ,lg+1) + + owx * wy * wz * Bzc(jg ,kg-1,lg+1) + + wx * wy * wz * Bzc(jg-1,kg-1,lg+1); + + // interp from fine stagged to fine nodal + Real bf = 0.25*(Bzf(j-1,k-1,l) + Bzf(j,k-1,l) + Bzf(j-1,k,l) + Bzf(j,k,l)); + +#endif + + Bza(j,k,l) = bg + (bf-bc); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bxa, + amrex::Array4<amrex::Real const> const& Bxf) +{ +#if (AMREX_SPACEDIM == 2) + Bxa(j,k,0) = 0.5*(Bxf(j,k-1,0) + Bxf(j,k,0)); +#else + Bxa(j,k,l) = 0.25*(Bxf(j,k-1,l-1) + Bxf(j,k,l-1) + Bxf(j,k-1,l) + Bxf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bya, + amrex::Array4<amrex::Real const> const& Byf) +{ +#if (AMREX_SPACEDIM == 2) + Bya(j,k,0) = 0.25*(Byf(j,k,0) + Byf(j-1,k,0) + Byf(j,k-1,0) + Byf(j-1,k-1,0)); +#else + Bya(j,k,l) = 0.25*(Byf(j-1,k,l-1) + Byf(j,k,l-1) + Byf(j-1,k,l) + Byf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bza, + amrex::Array4<amrex::Real const> const& Bzf) +{ +#if (AMREX_SPACEDIM == 2) + Bza(j,k,0) = 0.5*(Bzf(j-1,k,0) + Bzf(j,k,0)); +#else + Bza(j,k,l) = 0.25*(Bzf(j-1,k-1,l) + Bzf(j,k-1,l) + Bzf(j-1,k,l) + Bzf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Exa, + amrex::Array4<amrex::Real const> const& Exf, + amrex::Array4<amrex::Real const> const& Exc, + amrex::Array4<amrex::Real const> const& Exg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * Exg(jg ,kg ,0) + + owx * wy * Exg(jg ,kg+1,0) + + wx * owy * Exg(jg+1,kg ,0) + + wx * wy * Exg(jg+1,kg+1,0); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + Real ec = owx * owy * Exc(jg ,kg ,0) + + owx * wy * Exc(jg ,kg+1,0) + + wx * owy * Exc(jg-1,kg ,0) + + wx * wy * Exc(jg-1,kg+1,0); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Exf(j-1,k,0) + Exf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * owz * Exg(jg ,kg ,lg ) + + wx * owy * owz * Exg(jg+1,kg ,lg ) + + owx * wy * owz * Exg(jg ,kg+1,lg ) + + wx * wy * owz * Exg(jg+1,kg+1,lg ) + + owx * owy * wz * Exg(jg ,kg ,lg+1) + + wx * owy * wz * Exg(jg+1,kg ,lg+1) + + owx * wy * wz * Exg(jg ,kg+1,lg+1) + + wx * wy * wz * Exg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + Real ec = owx * owy * owz * Exc(jg ,kg ,lg ) + + wx * owy * owz * Exc(jg-1,kg ,lg ) + + owx * wy * owz * Exc(jg ,kg+1,lg ) + + wx * wy * owz * Exc(jg-1,kg+1,lg ) + + owx * owy * wz * Exc(jg ,kg ,lg+1) + + wx * owy * wz * Exc(jg-1,kg ,lg+1) + + owx * wy * wz * Exc(jg ,kg+1,lg+1) + + wx * wy * wz * Exc(jg-1,kg+1,lg+1); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Exf(j-1,k,l) + Exf(j,k,l)); + +#endif + + Exa(j,k,l) = eg + (ef-ec); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eya, + amrex::Array4<amrex::Real const> const& Eyf, + amrex::Array4<amrex::Real const> const& Eyc, + amrex::Array4<amrex::Real const> const& Eyg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal and coarse staggered to fine nodal + Real eg = owx * owy * (Eyg(jg ,kg ,0) + Eyc(jg ,kg ,0)) + + owx * wy * (Eyg(jg ,kg+1,0) + Eyc(jg ,kg+1,0)) + + wx * owy * (Eyg(jg+1,kg ,0) + Eyc(jg+1,kg ,0)) + + wx * wy * (Eyg(jg+1,kg+1,0) + Eyc(jg+1,kg+1,0)); + Real ec = 0.0; + + // interp from fine staggered to fine nodal + Real ef = Eyf(j,k,0); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * owz * Eyg(jg ,kg ,lg ) + + wx * owy * owz * Eyg(jg+1,kg ,lg ) + + owx * wy * owz * Eyg(jg ,kg+1,lg ) + + wx * wy * owz * Eyg(jg+1,kg+1,lg ) + + owx * owy * wz * Eyg(jg ,kg ,lg+1) + + wx * owy * wz * Eyg(jg+1,kg ,lg+1) + + owx * wy * wz * Eyg(jg ,kg+1,lg+1) + + wx * wy * wz * Eyg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + Real ec = owx * owy * owz * Eyc(jg ,kg ,lg ) + + wx * owy * owz * Eyc(jg+1,kg ,lg ) + + owx * wy * owz * Eyc(jg ,kg-1,lg ) + + wx * wy * owz * Eyc(jg+1,kg-1,lg ) + + owx * owy * wz * Eyc(jg ,kg ,lg+1) + + wx * owy * wz * Eyc(jg+1,kg ,lg+1) + + owx * wy * wz * Eyc(jg ,kg-1,lg+1) + + wx * wy * wz * Eyc(jg+1,kg-1,lg+1); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Eyf(j,k-1,l) + Eyf(j,k,l)); + +#endif + + Eya(j,k,l) = eg + (ef-ec); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eza, + amrex::Array4<amrex::Real const> const& Ezf, + amrex::Array4<amrex::Real const> const& Ezc, + amrex::Array4<amrex::Real const> const& Ezg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * Ezg(jg ,kg ,0) + + owx * wy * Ezg(jg ,kg+1,0) + + wx * owy * Ezg(jg+1,kg ,0) + + wx * wy * Ezg(jg+1,kg+1,0); + + // interp from coarse stagged to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + Real ec = owx * owy * Ezc(jg ,kg ,0) + + owx * wy * Ezc(jg ,kg-1,0) + + wx * owy * Ezc(jg+1,kg ,0) + + wx * wy * Ezc(jg+1,kg-1,0); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Ezf(j,k-1,0) + Ezf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * owz * Ezg(jg ,kg ,lg ) + + wx * owy * owz * Ezg(jg+1,kg ,lg ) + + owx * wy * owz * Ezg(jg ,kg+1,lg ) + + wx * wy * owz * Ezg(jg+1,kg+1,lg ) + + owx * owy * wz * Ezg(jg ,kg ,lg+1) + + wx * owy * wz * Ezg(jg+1,kg ,lg+1) + + owx * wy * wz * Ezg(jg ,kg+1,lg+1) + + wx * wy * wz * Ezg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wz = 0.5-wz; owz = 1.0-wz; + Real ec = owx * owy * owz * Ezc(jg ,kg ,lg ) + + wx * owy * owz * Ezc(jg+1,kg ,lg ) + + owx * wy * owz * Ezc(jg ,kg+1,lg ) + + wx * wy * owz * Ezc(jg+1,kg+1,lg ) + + owx * owy * wz * Ezc(jg ,kg ,lg-1) + + wx * owy * wz * Ezc(jg+1,kg ,lg-1) + + owx * wy * wz * Ezc(jg ,kg+1,lg-1) + + wx * wy * wz * Ezc(jg+1,kg+1,lg-1); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Ezf(j,k,l-1) + Ezf(j,k,l)); + +#endif + + Eza(j,k,l) = eg + (ef-ec); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Exa, + amrex::Array4<amrex::Real const> const& Exf) +{ + Exa(j,k,l) = 0.5*(Exf(j-1,k,l) + Exf(j,k,l)); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eya, + amrex::Array4<amrex::Real const> const& Eyf) +{ +#if (AMREX_SPACEDIM == 2) + Eya(j,k,0) = Eyf(j,k,0); +#else + Eya(j,k,l) = 0.5*(Eyf(j,k-1,l) + Eyf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eza, + amrex::Array4<amrex::Real const> const& Ezf) +{ +#if (AMREX_SPACEDIM == 2) + Eza(j,k,0) = 0.5*(Ezf(j,k-1,0) + Ezf(j,k,0)); +#else + Eza(j,k,l) = 0.5*(Ezf(j,k,l-1) + Ezf(j,k,l)); +#endif +} + #endif diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index 5441755f5..2ae167283 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -91,7 +91,7 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa // Aux patch - if (lev == 0) + if (lev == 0 && Bfield_aux[0][0]->ixType() == Bfield_fp[0][0]->ixType()) { for (int idim = 0; idim < 3; ++idim) { Bfield_aux[lev][idim].reset(new MultiFab(*Bfield_fp[lev][idim], amrex::make_alias, 0, Bfield_aux[lev][idim]->nComp())); |