aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/BoundaryConditions/WarpX_PML_kernels.H2
-rw-r--r--Source/Diagnostics/FieldIO.H2
-rw-r--r--Source/Parallelization/Make.package1
-rw-r--r--Source/Parallelization/WarpXComm.cpp147
-rw-r--r--Source/Parallelization/WarpXComm_K.H161
-rw-r--r--Source/Parallelization/WarpXSumGuardCells.H6
6 files changed, 218 insertions, 101 deletions
diff --git a/Source/BoundaryConditions/WarpX_PML_kernels.H b/Source/BoundaryConditions/WarpX_PML_kernels.H
index 23d19e2e8..8f779a5c2 100644
--- a/Source/BoundaryConditions/WarpX_PML_kernels.H
+++ b/Source/BoundaryConditions/WarpX_PML_kernels.H
@@ -99,7 +99,7 @@ void warpx_push_pml_bx_ckc(int i, int j, int k, Array4<Real> const&Bx,
- Ey(i+1,j+1,k ,0) - Ey(i+1,j+1,k ,1) - Ey(i+1,j+1,k ,2)
+ Ey(i-1,j+1,k+1,0) + Ey(i-1,j+1,k+1,1) + Ey(i-1,j+1,k+1,2)
- Ey(i-1,j+1,k ,0) - Ey(i-1,j+1,k ,1) - Ey(i-1,j+1,k ,2)
- + Ey(i+1,j-1,k+1,0) + Ey(i+1,j-1,k+1,1) + Ey(i+1,j-2,k+1,2)
+ + Ey(i+1,j-1,k+1,0) + Ey(i+1,j-1,k+1,1) + Ey(i+1,j-1,k+1,2)
- Ey(i+1,j-1,k ,0) - Ey(i+1,j-1,k ,1) - Ey(i+1,j-1,k ,2)
+ Ey(i-1,j-1,k+1,0) + Ey(i-1,j-1,k+1,1) + Ey(i-1,j-1,k+1,1)
- Ey(i-1,j-1,k ,0) - Ey(i-1,j-1,k ,1) - Ey(i-1,j-1,k ,2));
diff --git a/Source/Diagnostics/FieldIO.H b/Source/Diagnostics/FieldIO.H
index f08d85f2d..7cdc9b710 100644
--- a/Source/Diagnostics/FieldIO.H
+++ b/Source/Diagnostics/FieldIO.H
@@ -99,7 +99,7 @@ getReversedVec( const amrex::Real* v );
void
WriteOpenPMDFields( const std::string& filename,
const std::vector<std::string>& varnames,
- const amrex::MultiFab& mf, const Geometry& geom,
+ const amrex::MultiFab& mf, const amrex::Geometry& geom,
const int iteration, const double time );
#endif // WARPX_USE_OPENPMD
diff --git a/Source/Parallelization/Make.package b/Source/Parallelization/Make.package
index 3d1fcf1da..c74583522 100644
--- a/Source/Parallelization/Make.package
+++ b/Source/Parallelization/Make.package
@@ -1,6 +1,7 @@
CEXE_sources += WarpXComm.cpp
CEXE_sources += WarpXRegrid.cpp
CEXE_headers += WarpXSumGuardCells.H
+CEXE_headers += WarpXComm_K.H
INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Parallelization
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Parallelization
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index e24dd772c..990d0f988 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -1,9 +1,8 @@
+#include <WarpXComm_K.H>
#include <WarpX.H>
#include <WarpX_f.H>
#include <WarpXSumGuardCells.H>
-#include <AMReX_FillPatchUtil_F.H>
-
#include <algorithm>
#include <cstdlib>
@@ -52,8 +51,6 @@ WarpX::UpdateAuxilaryData ()
{
BL_PROFILE("UpdateAuxilaryData()");
- const int use_limiter = 0;
-
for (int lev = 1; lev <= finest_level; ++lev)
{
const auto& crse_period = Geom(lev-1).periodicity();
@@ -81,57 +78,37 @@ WarpX::UpdateAuxilaryData ()
MultiFab::Subtract(dBy, *Bfield_cp[lev][1], 0, 0, Bfield_cp[lev][1]->nComp(), ng);
MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, Bfield_cp[lev][2]->nComp(), ng);
- const Real* dx = Geom(lev-1).CellSize();
const int refinement_ratio = refRatio(lev-1)[0];
+ AMREX_ALWAYS_ASSERT(refinement_ratio == 2);
+
#ifdef _OPENMP
-#pragma omp parallel
+#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
+ for (MFIter mfi(*Bfield_aux[lev][0]); mfi.isValid(); ++mfi)
{
- std::array<FArrayBox,3> bfab;
- 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_c = dBx.const_array(mfi);
+ Array4<Real const> const& by_c = dBy.const_array(mfi);
+ Array4<Real const> const& bz_c = dBz.const_array(mfi);
+
+ amrex::ParallelFor(Box(bx_aux), Box(by_aux), Box(bz_aux),
+ [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
{
- Box ccbx = mfi.fabbox();
- ccbx.enclosedCells();
- ccbx.coarsen(refinement_ratio).refine(refinement_ratio); // so that ccbx is coarsenable
-
- const FArrayBox& cxfab = dBx[mfi];
- const FArrayBox& cyfab = dBy[mfi];
- const FArrayBox& czfab = dBz[mfi];
- bfab[0].resize(amrex::convert(ccbx,Bx_nodal_flag));
- bfab[1].resize(amrex::convert(ccbx,By_nodal_flag));
- bfab[2].resize(amrex::convert(ccbx,Bz_nodal_flag));
-
-#if (AMREX_SPACEDIM == 3)
- amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(),
- BL_TO_FORTRAN_ANYD(bfab[0]),
- BL_TO_FORTRAN_ANYD(bfab[1]),
- BL_TO_FORTRAN_ANYD(bfab[2]),
- BL_TO_FORTRAN_ANYD(cxfab),
- BL_TO_FORTRAN_ANYD(cyfab),
- BL_TO_FORTRAN_ANYD(czfab),
- dx, &refinement_ratio,&use_limiter);
-#else
- amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(),
- BL_TO_FORTRAN_ANYD(bfab[0]),
- BL_TO_FORTRAN_ANYD(bfab[2]),
- BL_TO_FORTRAN_ANYD(cxfab),
- BL_TO_FORTRAN_ANYD(czfab),
- dx, &refinement_ratio,&use_limiter);
- amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(),
- BL_TO_FORTRAN_ANYD(bfab[1]),
- BL_TO_FORTRAN_ANYD(cyfab),
- &refinement_ratio,&use_limiter);
-#endif
-
- for (int idim = 0; idim < 3; ++idim)
- {
- FArrayBox& aux = (*Bfield_aux[lev][idim])[mfi];
- FArrayBox& fp = (*Bfield_fp[lev][idim])[mfi];
- const Box& bx = aux.box();
- aux.copy(fp, bx, 0, bx, 0, 1);
- aux.plus(bfab[idim], bx, bx, 0, 0, 1);
- }
- }
+ warpx_interp_bfield_x(j,k,l, bx_aux, bx_fp, bx_c);
+ },
+ [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+ {
+ warpx_interp_bfield_y(j,k,l, by_aux, by_fp, by_c);
+ },
+ [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+ {
+ warpx_interp_bfield_z(j,k,l, bz_aux, bz_fp, bz_c);
+ });
}
}
@@ -156,56 +133,34 @@ WarpX::UpdateAuxilaryData ()
MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, Efield_cp[lev][1]->nComp(), ng);
MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, Efield_cp[lev][2]->nComp(), ng);
- const int refinement_ratio = refRatio(lev-1)[0];
#ifdef _OPEMP
-#pragma omp parallel
+#pragma omp parallel if (Gpu::notInLaunchRegion())
#endif
+ for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi)
{
- std::array<FArrayBox,3> efab;
- 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_c = dEx.const_array(mfi);
+ Array4<Real const> const& ey_c = dEy.const_array(mfi);
+ Array4<Real const> const& ez_c = dEz.const_array(mfi);
+
+ amrex::ParallelFor(Box(ex_aux), Box(ey_aux), Box(ez_aux),
+ [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
{
- Box ccbx = mfi.fabbox();
- ccbx.enclosedCells();
- ccbx.coarsen(refinement_ratio).refine(refinement_ratio); // so that ccbx is coarsenable
-
- const FArrayBox& cxfab = dEx[mfi];
- const FArrayBox& cyfab = dEy[mfi];
- const FArrayBox& czfab = dEz[mfi];
- efab[0].resize(amrex::convert(ccbx,Ex_nodal_flag));
- efab[1].resize(amrex::convert(ccbx,Ey_nodal_flag));
- efab[2].resize(amrex::convert(ccbx,Ez_nodal_flag));
-
-#if (AMREX_SPACEDIM == 3)
- amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(),
- BL_TO_FORTRAN_ANYD(efab[0]),
- BL_TO_FORTRAN_ANYD(efab[1]),
- BL_TO_FORTRAN_ANYD(efab[2]),
- BL_TO_FORTRAN_ANYD(cxfab),
- BL_TO_FORTRAN_ANYD(cyfab),
- BL_TO_FORTRAN_ANYD(czfab),
- &refinement_ratio,&use_limiter);
-#else
- amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(),
- BL_TO_FORTRAN_ANYD(efab[0]),
- BL_TO_FORTRAN_ANYD(efab[2]),
- BL_TO_FORTRAN_ANYD(cxfab),
- BL_TO_FORTRAN_ANYD(czfab),
- &refinement_ratio,&use_limiter);
- amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(),
- BL_TO_FORTRAN_ANYD(efab[1]),
- BL_TO_FORTRAN_ANYD(cyfab),
- &refinement_ratio);
-#endif
-
- for (int idim = 0; idim < 3; ++idim)
- {
- FArrayBox& aux = (*Efield_aux[lev][idim])[mfi];
- FArrayBox& fp = (*Efield_fp[lev][idim])[mfi];
- const Box& bx = aux.box();
- aux.copy(fp, bx, 0, bx, 0, Efield_fp[lev][idim]->nComp());
- aux.plus(efab[idim], bx, bx, 0, 0, Efield_fp[lev][idim]->nComp());
- }
- }
+ warpx_interp_efield_x(j,k,l, ex_aux, ex_fp, ex_c);
+ },
+ [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+ {
+ warpx_interp_efield_y(j,k,l, ey_aux, ey_fp, ey_c);
+ },
+ [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept
+ {
+ warpx_interp_efield_z(j,k,l, ez_aux, ez_fp, ez_c);
+ });
}
}
}
diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H
new file mode 100644
index 000000000..093323ec3
--- /dev/null
+++ b/Source/Parallelization/WarpXComm_K.H
@@ -0,0 +1,161 @@
+#ifndef WARPX_COMM_K_H_
+#define WARPX_COMM_K_H_
+
+#include <AMReX_FArrayBox.H>
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_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)
+{
+ using namespace amrex;
+
+ int lg = amrex::coarsen(l,2);
+ int kg = amrex::coarsen(k,2);
+ int jg = amrex::coarsen(j,2);
+
+ Real wx = (j == jg*2) ? 0.0 : 0.5;
+ Real owx = 1.0-wx;
+ Bxa(j,k,l) = owx * Bxc(jg,kg,lg) + wx * Bxc(jg+1,kg,lg) + Bxf(j,k,l);
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_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)
+{
+ using namespace amrex;
+
+ int lg = amrex::coarsen(l,2);
+ int kg = amrex::coarsen(k,2);
+ int jg = amrex::coarsen(j,2);
+
+ // Note that for 2d, l=0, because the amrex convention is used here.
+
+#if (AMREX_SPACEDIM == 3)
+ Real wy = (k == kg*2) ? 0.0 : 0.5;
+ Real owy = 1.0-wy;
+ Bya(j,k,l) = owy * Byc(jg,kg,lg) + wy * Byc(jg,kg+1,lg) + Byf(j,k,l);
+#else
+ Bya(j,k,l) = Byc(jg,kg,lg) + Byf(j,k,l);
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_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)
+{
+ using namespace amrex;
+
+ int lg = amrex::coarsen(l,2);
+ int kg = amrex::coarsen(k,2);
+ int jg = amrex::coarsen(j,2);
+
+ // Note that for 2d, l=0, because the amrex convention is used here.
+
+#if (AMREX_SPACEDIM == 3)
+ Real wz = (l == lg*2) ? 0.0 : 0.5;
+ Real owz = 1.0-wz;
+ Bza(j,k,l) = owz * Bzc(jg,kg,lg) + owz * Bzc(jg,kg,lg+1) + Bzf(j,k,l);
+#else
+ Real wy = (k == kg*2) ? 0.0 : 0.5;
+ Real owy = 1.0-wy;
+ Bza(j,k,l) = owy * Bzc(jg,kg,lg) + owy * Bzc(jg,kg+1,lg) + Bzf(j,k,l);
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_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)
+{
+ using namespace amrex;
+
+ int lg = amrex::coarsen(l,2);
+ int kg = amrex::coarsen(k,2);
+ int jg = amrex::coarsen(j,2);
+
+ Real wy = (k == kg*2) ? 0.0 : 0.5;
+ Real owy = 1.0-wy;
+
+#if (AMREX_SPACEDIM == 3)
+ Real wz = (l == lg*2) ? 0.0 : 0.5;
+ Real owz = 1.0-wz;
+ Exa(j,k,l) = owy * owz * Exc(jg ,kg ,lg )
+ + wy * owz * Exc(jg ,kg+1,lg )
+ + owy * wz * Exc(jg ,kg ,lg+1)
+ + wy * wz * Exc(jg ,kg+1,lg+1)
+ + Exf(j,k,l);
+#else
+ Exa(j,k,l) = owy * Exc(jg,kg,lg) + wy * Exc(jg,kg+1,lg) + Exf(j,k,l);
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_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)
+{
+ using namespace amrex;
+
+ int lg = amrex::coarsen(l,2);
+ int kg = amrex::coarsen(k,2);
+ int jg = amrex::coarsen(j,2);
+
+ Real wx = (j == jg*2) ? 0.0 : 0.5;
+ Real owx = 1.0-wx;
+
+#if (AMREX_SPACEDIM == 3)
+ Real wz = (l == lg*2) ? 0.0 : 0.5;
+ Real owz = 1.0-wz;
+ Eya(j,k,l) = owx * owz * Eyc(jg ,kg ,lg )
+ + wx * owz * Eyc(jg+1,kg ,lg )
+ + owx * wz * Eyc(jg ,kg ,lg+1)
+ + wx * wz * Eyc(jg+1,kg ,lg+1)
+ + Eyf(j,k,l);
+#else
+ Real wy = (k == kg*2) ? 0.0 : 0.5;
+ Real owy = 1.0-wy;
+ Eya(j,k,l) = owx * owy * Eyc(jg ,kg ,lg)
+ + wx * owy * Eyc(jg+1,kg ,lg)
+ + owx * wy * Eyc(jg ,kg+1,lg)
+ + wx * wy * Eyc(jg+1,kg+1,lg)
+ + Eyf(j,k,l);
+#endif
+}
+
+AMREX_GPU_DEVICE AMREX_FORCE_INLINE
+void warpx_interp_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)
+{
+ using namespace amrex;
+
+ int lg = amrex::coarsen(l,2);
+ int kg = amrex::coarsen(k,2);
+ int jg = amrex::coarsen(j,2);
+
+ Real wx = (j == jg*2) ? 0.0 : 0.5;
+ Real owx = 1.0-wx;
+
+#if (AMREX_SPACEDIM == 3)
+ Real wy = (k == kg*2) ? 0.0 : 0.5;
+ Real owy = 1.0-wy;
+ Eza(j,k,l) = owx * owy * Ezc(jg ,kg ,lg )
+ + wx * owy * Ezc(jg+1,kg ,lg )
+ + owx * wy * Ezc(jg ,kg+1,lg )
+ + wx * wy * Ezc(jg+1,kg+1,lg )
+ + Ezf(j,k,l);
+#else
+ Eza(j,k,l) = owx * Ezc(jg,kg,lg) + wx * Ezc(jg+1,kg,lg) + Ezf(j,k,l);
+#endif
+}
+
+#endif
diff --git a/Source/Parallelization/WarpXSumGuardCells.H b/Source/Parallelization/WarpXSumGuardCells.H
index 24ad1b80f..ce353c2b6 100644
--- a/Source/Parallelization/WarpXSumGuardCells.H
+++ b/Source/Parallelization/WarpXSumGuardCells.H
@@ -15,7 +15,7 @@
* updates both the *valid* cells and *guard* cells. (This is because a
* spectral solver requires the value of the sources over a large stencil.)
*/
-void
+inline void
WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
const int icomp=0, const int ncomp=1){
#ifdef WARPX_USE_PSATD
@@ -43,7 +43,7 @@ WarpXSumGuardCells(amrex::MultiFab& mf, const amrex::Periodicity& period,
* Note: `i_comp` is the component where the results will be stored in `dst`;
* The component from which we copy in `src` is always 0.
*/
-void
+inline void
WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src,
const amrex::Periodicity& period,
const int icomp=0, const int ncomp=1){
@@ -54,7 +54,7 @@ WarpXSumGuardCells(amrex::MultiFab& dst, amrex::MultiFab& src,
// Update only the valid cells
const amrex::IntVect n_updated_guards = amrex::IntVect::TheZeroVector();
#endif
- src.SumBoundary(icomp, ncomp, n_updated_guards, period);
+ src.SumBoundary(0, ncomp, n_updated_guards, period);
amrex::Copy( dst, src, 0, icomp, ncomp, n_updated_guards );
}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-17 10:56:18 +0000