Multilevel version of ParticleToMesh (AMReX-Codes#1945)

This uses the same algorithm as the AssignDensity method in `AmrCore/AMReX_AmrParticles.H` to handle the coarse / fine boundaries. However, it allows a user-specified lambda function to define the interpolation operation, as in the single-level `ParticleToMesh` routine.
cdaley · Apr 15, 2021 · d2e37b0 · d2e37b0
1 parent b1d05de
commit d2e37b0
Show file tree

Hide file tree

Showing 9 changed files with 384 additions and 9 deletions.
diff --git a/Src/AmrCore/AMReX_AmrParticles.H b/Src/AmrCore/AMReX_AmrParticles.H
@@ -14,9 +14,9 @@ template <int NStructReal, int NStructInt, int NArrayReal, int NArrayInt,
           template<class> class Allocator>
 void
 ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>
-::AssignDensity(int rho_index,
-                Vector<std::unique_ptr<MultiFab> >& mf_to_be_filled,
-                int lev_min, int ncomp, int finest_level, int ngrow) const
+::AssignDensity (int rho_index,
+                 Vector<std::unique_ptr<MultiFab> >& mf_to_be_filled,
+                 int lev_min, int ncomp, int finest_level, int ngrow) const
 {
 
     BL_PROFILE("ParticleContainer::AssignDensity()");
@@ -151,6 +151,100 @@ ParticleContainer<NStructReal, NStructInt, NArrayReal, NArrayInt, Allocator>
     }
 }
 
+template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
+void
+ParticleToMesh (PC const& pc, const Vector<MultiFab*>& mf,
+                int lev_min, int lev_max, F&& f,
+                bool zero_out_input=true, bool vol_weight=true)
+{
+    BL_PROFILE("amrex::ParticleToMesh");
+
+    if (lev_max == -1) { lev_max = pc.finestLevel(); }
+    while (!pc.GetParGDB()->LevelDefined(lev_max)) { lev_max--; }
+
+    if (lev_max == 0)
+    {
+        ParticleToMesh(pc, *mf[0], 0, std::forward<F>(f), zero_out_input);
+        if (vol_weight) {
+            const Real* dx = pc.Geom(0).CellSize();
+            const Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]);
+            mf[0]->mult(1.0/vol, 0, mf[0]->nComp(), mf[0]->nGrow());
+        }
+        return;
+    }
+
+    if (zero_out_input)
+    {
+        for (int lev = lev_min; lev <= lev_max; ++lev)
+        {
+            mf[lev]->setVal(0.0);
+        }
+    }
+
+    Vector<MultiFab> mf_part(lev_max+1);
+    Vector<MultiFab> mf_tmp( lev_max+1);
+    for (int lev = lev_min; lev <= lev_max; ++lev)
+    {
+        mf_part[lev].define(pc.ParticleBoxArray(lev),
+                           pc.ParticleDistributionMap(lev),
+                           mf[lev]->nComp(), mf[lev]->nGrowVect());
+        mf_tmp[lev].define( pc.ParticleBoxArray(lev),
+                            pc.ParticleDistributionMap(lev),
+                            mf[lev]->nComp(), mf[lev]->nGrowVect());
+        mf_part[lev].setVal(0.0);
+        mf_tmp[lev].setVal( 0.0);
+    }
+
+    // configure this to do a no-op at the physical boundaries.
+    int lo_bc[] = {BCType::int_dir, BCType::int_dir, BCType::int_dir};
+    int hi_bc[] = {BCType::int_dir, BCType::int_dir, BCType::int_dir};
+    Vector<BCRec> bcs(1, BCRec(lo_bc, hi_bc));
+    PCInterp mapper;
+
+    for (int lev = lev_min; lev <= lev_max; ++lev)
+    {
+        ParticleToMesh(pc, mf_part[lev], lev, std::forward<F>(f), zero_out_input);
+        if (vol_weight) {
+            const Real* dx = pc.Geom(lev).CellSize();
+            const Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]);
+            mf_part[lev].mult(1.0/vol, 0, mf_part[lev].nComp(), mf_part[lev].nGrow());
+        }
+
+        if (lev < lev_max) {
+            PhysBCFunctNoOp cphysbc, fphysbc;
+            amrex::InterpFromCoarseLevel(mf_tmp[lev+1], 0.0, mf_part[lev],
+                                         0, 0, mf_part[lev].nComp(),
+                                         pc.GetParGDB()->Geom(lev),
+                                         pc.GetParGDB()->Geom(lev+1),
+                                         cphysbc, 0, fphysbc, 0,
+                                         pc.GetParGDB()->refRatio(lev), &mapper,
+                                         bcs, 0);
+        }
+
+        if (lev > lev_min) {
+            // Note - this will double count the mass on the coarse level in
+            // regions covered by the fine level, but this will be corrected
+            // below in the call to average_down.
+            amrex::sum_fine_to_coarse(mf_part[lev], mf_part[lev-1], 0, mf_part[lev].nComp(),
+                                      pc.GetParGDB()->refRatio(lev-1),
+                                      pc.GetParGDB()->Geom(lev-1),
+                                      pc.GetParGDB()->Geom(lev));
+        }
+
+        mf_part[lev].plus(mf_tmp[lev], 0, mf_part[lev].nComp(), 0);
+    }
+
+    for (int lev = lev_max-1; lev >= lev_min; --lev) {
+        amrex::average_down(mf_part[lev+1], mf_part[lev], 0, mf_part[lev].nComp(),
+                            pc.GetParGDB()->refRatio(lev));
+    }
+
+    for (int lev = lev_min; lev <= lev_max; lev++)
+    {
+        mf[lev]->ParallelCopy(mf_part[lev],0,0,mf_part[lev].nComp(),0,0);
+    }
+}
+
 template <int NStructReal, int NStructInt=0, int NArrayReal=0, int NArrayInt=0,
           template<class> class Allocator=DefaultAllocator>
 class AmrParticleContainer

diff --git a/Src/Particle/AMReX_ParticleMesh.H b/Src/Particle/AMReX_ParticleMesh.H
@@ -4,6 +4,7 @@
 
 #include <AMReX_TypeTraits.H>
 #include <AMReX_MultiFab.H>
+#include <AMReX_ParticleUtil.H>
 
 namespace amrex
 {
@@ -21,8 +22,10 @@ ParticleToMesh (PC const& pc, MF& mf, int lev, F&& f, bool zero_out_input=true)
                     mf.nComp(), mf.nGrowVect());
     mf_tmp.setVal(0.0);
 
+    const auto plo = pc.Geom(lev).ProbLoArray();
+    const auto dxi = pc.Geom(lev).InvCellSizeArray();
+
     using ParIter = typename PC::ParConstIterType;
-    const auto& plevel = pc.GetParticles(lev);
 #ifdef AMREX_USE_GPU
     if (Gpu::inLaunchRegion())
     {
@@ -38,7 +41,7 @@ ParticleToMesh (PC const& pc, MF& mf, int lev, F&& f, bool zero_out_input=true)
 
             AMREX_FOR_1D( np, i,
             {
-                f(pstruct[i], fabarr);
+                particle_detail::call_f(f, pstruct[i], fabarr, plo, dxi);
             });
         }
     }
@@ -67,7 +70,7 @@ ParticleToMesh (PC const& pc, MF& mf, int lev, F&& f, bool zero_out_input=true)
 
                 AMREX_FOR_1D( np, i,
                 {
-                    f(pstruct[i], fabarr);
+                    particle_detail::call_f(f, pstruct[i], fabarr, plo, dxi);
                 });
 
                 fab.atomicAdd<RunOn::Host>(local_fab, tile_box, tile_box, 0, 0, mf_tmp.nComp());
@@ -91,7 +94,9 @@ MeshToParticle (PC& pc, MF const& mf, int lev, F&& f)
 
     if (mf_pointer != &mf) mf_pointer->copy(mf,0,0,mf.nComp(),mf.nGrowVect(),mf.nGrowVect());
 
-    auto& plevel = pc.GetParticles(lev);
+    const auto plo = pc.Geom(lev).ProbLoArray();
+    const auto dxi = pc.Geom(lev).InvCellSizeArray();
+
     using ParIter = typename PC::ParIterType;
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
@@ -108,7 +113,7 @@ MeshToParticle (PC& pc, MF const& mf, int lev, F&& f)
 
         AMREX_FOR_1D( np, i,
         {
-            f(pstruct[i], fabarr);
+            particle_detail::call_f(f, pstruct[i], fabarr, plo, dxi);
         });
     }
 

diff --git a/Src/Particle/AMReX_ParticleUtil.H b/Src/Particle/AMReX_ParticleUtil.H
@@ -53,6 +53,46 @@ auto call_f (F const& f, SrcData const& src, N i, amrex::RandomEngine const&) no
     return f(src,i);
 }
 
+template <typename F, typename P, typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+auto call_f (F const& f, P const& p, Array4<T> const& fabarr,
+             GpuArray<Real,AMREX_SPACEDIM> const& plo,
+             GpuArray<Real,AMREX_SPACEDIM> const& dxi) noexcept
+    -> decltype(f(p, fabarr, plo, dxi))
+{
+    return f(p, fabarr, plo, dxi);
+}
+
+template <typename F, typename P, typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+auto call_f (F const& f, P const& p, Array4<T> const& fabarr,
+             GpuArray<Real,AMREX_SPACEDIM> const&,
+             GpuArray<Real,AMREX_SPACEDIM> const&) noexcept
+    -> decltype(f(p, fabarr))
+{
+    return f(p, fabarr);
+}
+
+template <typename F, typename P, typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+auto call_f (F const& f, P& p, Array4<const T> const& fabarr,
+             GpuArray<Real,AMREX_SPACEDIM> const& plo,
+             GpuArray<Real,AMREX_SPACEDIM> const& dxi) noexcept
+    -> decltype(f(p, fabarr, plo, dxi))
+{
+    return f(p, fabarr, plo, dxi);
+}
+
+template <typename F, typename P, typename T>
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+auto call_f (F const& f, P& p, Array4<const T> const& fabarr,
+             GpuArray<Real,AMREX_SPACEDIM> const&,
+             GpuArray<Real,AMREX_SPACEDIM> const&) noexcept
+    -> decltype(f(p, fabarr))
+{
+    return f(p, fabarr);
+}
+
 }
 
 /**

diff --git a/Tests/Particles/ParticleMesh/main.cpp b/Tests/Particles/ParticleMesh/main.cpp
@@ -18,7 +18,7 @@ struct TestParams {
   bool verbose;
 };
 
-void testParticleMesh(TestParams& parms)
+void testParticleMesh (TestParams& parms)
 {
 
   RealBox real_box;

diff --git a/Tests/Particles/ParticleMeshMultiLevel/CMakeLists.txt b/Tests/Particles/ParticleMeshMultiLevel/CMakeLists.txt
@@ -0,0 +1,7 @@
+set(_sources     main.cpp)
+set(_input_files inputs  )
+
+setup_test(_sources _input_files)
+
+unset(_sources)
+unset(_input_files)
diff --git a/Tests/Particles/ParticleMeshMultiLevel/GNUmakefile b/Tests/Particles/ParticleMeshMultiLevel/GNUmakefile
@@ -0,0 +1,30 @@
+AMREX_HOME = ../../../
+
+DEBUG	= TRUE
+DEBUG	= FALSE
+
+DIM	= 3
+
+COMP    = gcc
+
+TINY_PROFILE = TRUE
+USE_PARTICLES = TRUE
+
+PRECISION = DOUBLE
+
+USE_MPI   = TRUE
+USE_OMP   = FALSE
+
+###################################################
+
+EBASE     = main
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.defs
+
+include ./Make.package
+include $(AMREX_HOME)/Src/Base/Make.package
+include $(AMREX_HOME)/Src/Boundary/Make.package
+include $(AMREX_HOME)/Src/Particle/Make.package
+include $(AMREX_HOME)/Src/AmrCore/Make.package
+
+include $(AMREX_HOME)/Tools/GNUMake/Make.rules
diff --git a/Tests/Particles/ParticleMeshMultiLevel/Make.package b/Tests/Particles/ParticleMeshMultiLevel/Make.package
@@ -0,0 +1,2 @@
+CEXE_sources += main.cpp
+
diff --git a/Tests/Particles/ParticleMeshMultiLevel/inputs b/Tests/Particles/ParticleMeshMultiLevel/inputs
@@ -0,0 +1,26 @@
+
+# Domain size
+
+#nx = 32 # number of grid points along the x axis
+#ny = 32 # number of grid points along the y axis
+#nz = 32 # number of grid points along the z axis
+
+#nx = 64 # number of grid points along the x axis
+#ny = 64 # number of grid points along the y axis
+#nz = 64 # number of grid points along the z axis
+
+nx = 128 # number of grid points along the x axis
+ny = 128 # number of grid points along the y axis
+nz = 128 # number of grid points along the z axis
+
+# Maximum allowable size of each subdomain in the problem domain;
+#    this is used to decompose the domain for parallel calculations.
+max_grid_size = 32
+
+# Number of particles per cell
+nppc = 10
+
+# Verbosity
+verbose = true   # set to true to get more verbosity
+
+nlevs=2