Add particles (AMReX-Fluids#99)

src/Make.incflo

@@ -15,6 +15,11 @@ Bdirs += src/rheology
 Bdirs += src/setup
 Bdirs += src/utilities
 
+ifeq ($(USE_PARTICLES), TRUE)
+  DEFINES += -DINCFLO_USE_PARTICLES
+  Bdirs += src/particles
+endif
+
 ifeq ($(USE_EB), TRUE)
   Bdirs += src/embedded_boundaries
   USERSuffix += .EB
@@ -53,6 +58,10 @@ VPATH_LOCATIONS   += $(Blocs)
 #These are the directories in AMReX
 Pdirs := Base AmrCore Boundary LinearSolvers/MLMG
 
+ifeq ($(USE_PARTICLES), TRUE)
+  Pdirs += Particle
+endif
+
 ifeq ($(USE_EB), TRUE)
   Pdirs += EB
 endif

src/incflo.H

@@ -14,6 +14,10 @@
 #include <AMReX_EBMultiFabUtil.H>
 #endif
 
+#ifdef INCFLO_USE_PARTICLES
+#include "ParticleData.H"
+#endif
+
 #include <DiffusionTensorOp.H>
 #include <DiffusionScalarOp.H>
 
@@ -353,6 +357,7 @@ private:
     bool use_tensor_correction = false;
 
     // AMR / refinement settings
+    int m_refine_particles = 1;
     int m_refine_cutcells = 1;
     int m_regrid_int = -1;
 
@@ -541,6 +546,32 @@ private:
     int m_plt_error_p     = 0;
     int m_plt_error_mac_p = 0;
 
+#ifdef INCFLO_USE_PARTICLES
+    int m_plt_particle_count = 0;
+
+    // Particle container with all particle species
+    ParticleData particleData;
+
+    // variables and functions for tracers particles
+    bool  m_use_tracer_particles; /*!< tracer particles that advect with flow */
+
+    /*! Read tracer particles parameters */
+    void readTracerParticlesParams ();
+
+    /*! Initialize tracer particles */
+    void initializeTracerParticles (amrex::ParGDBBase* gdb
+#ifdef AMREX_USE_EB
+                            ,amrex::EBFArrayBoxFactory const& ebfact
+#endif
+                           );
+
+    /*! Evolve tracers and hydro particles */
+    void evolveTracerParticles (AMREX_D_DECL(amrex::Vector<amrex::MultiFab const*> const& u_mac,
+                                             amrex::Vector<amrex::MultiFab const*> const& v_mac,
+                                             amrex::Vector<amrex::MultiFab const*> const& w_mac));
+
+#endif
+
     struct LevelData {
         LevelData () = default;
         LevelData (amrex::BoxArray const& ba,

src/incflo.cpp

@@ -55,6 +55,20 @@ void incflo::InitData ()
         // with MakeNewLevelFromScratch.
         InitFromScratch(m_cur_time);
 
+#ifdef AMREX_USE_EB
+#ifdef INCFLO_USE_PARTICLES
+        const auto& ebfact = EBFactory(0);
+#endif
+#endif
+
+#ifdef INCFLO_USE_PARTICLES
+        initializeTracerParticles( (ParGDBBase*)GetParGDB()
+#ifdef AMREX_USE_EB
+                                  ,ebfact
+#endif
+                                 );
+#endif
+
 #ifdef AMREX_USE_EB
         InitialRedistribution();
 #endif
@@ -89,6 +103,10 @@ void incflo::InitData ()
         // Read starting configuration from chk file.
         ReadCheckpointFile();
 
+#ifdef INCFLO_USE_PARTICLES
+        particleData.Redistribute();
+#endif
+
         if (m_plotfile_on_restart)
         {
             WritePlotFile();

src/incflo_advance.cpp

@@ -66,6 +66,10 @@ void incflo::Advance()
         ApplyCorrector();
     }
 
+#ifdef INCFLO_USE_PARTICLES
+    particleData.Redistribute();
+#endif
+
 #if 0
     // This sums over all levels
     if (m_test_tracer_conservation) {

src/incflo_apply_corrector.cpp

@@ -468,6 +468,14 @@ void incflo::ApplyCorrector()
     bool incremental_projection = false;
     ApplyProjection(GetVecOfConstPtrs(density_nph), new_time,m_dt,incremental_projection);
 
+#ifdef INCFLO_USE_PARTICLES
+    // **************************************************************************************
+    // Update the particle positions
+    // **************************************************************************************
+    evolveTracerParticles(AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac),
+                                       GetVecOfConstPtrs(w_mac)));
+#endif
+
 #ifdef AMREX_USE_EB
     // **********************************************************************************************
     //

src/incflo_apply_predictor.cpp

@@ -444,6 +444,16 @@ void incflo::ApplyPredictor (bool incremental_projection)
     // **********************************************************************************************
     ApplyProjection(GetVecOfConstPtrs(density_nph),new_time,m_dt,incremental_projection);
 
+#ifdef INCFLO_USE_PARTICLES
+    // **************************************************************************************
+    // Update the particle positions
+    // **************************************************************************************
+    if (m_advection_type != "MOL") {
+        evolveTracerParticles(AMREX_D_DECL(GetVecOfConstPtrs(u_mac), GetVecOfConstPtrs(v_mac),
+                                           GetVecOfConstPtrs(w_mac)));
+    }
+#endif
+
 #ifdef AMREX_USE_EB
     // **********************************************************************************************
     //

src/incflo_regrid.cpp

@@ -108,6 +108,10 @@ void incflo::RemakeLevel (int lev, Real time, const BoxArray& ba,
 #else
     macproj = std::make_unique<Hydro::MacProjector>(Geom(0,finest_level));
 #endif
+
+#ifdef INCFLO_USE_PARTICLES
+    particleData.Redistribute();
+#endif
 }
 
 // Delete level data

src/incflo_tagging.cpp

@@ -8,7 +8,7 @@ using namespace amrex;
 
 // tag cells for refinement
 // overrides the pure virtual function in AmrCore
-void incflo::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/)
+void incflo::ErrorEst (int levc, TagBoxArray& tags, Real time, int /*ngrow*/)
 {
     BL_PROFILE("incflo::ErrorEst()");
 
@@ -45,30 +45,30 @@ void incflo::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/)
 
     const auto   tagval = TagBox::SET;
 
-    bool tag_rho = lev < rhoerr_v.size();
-    bool tag_gradrho = lev < gradrhoerr_v.size();
+    bool tag_rho = levc < rhoerr_v.size();
+    bool tag_gradrho = levc < gradrhoerr_v.size();
 
     if (tag_gradrho) {
-        fillpatch_density(lev, time, m_leveldata[lev]->density, 1);
+        fillpatch_density(levc, time, m_leveldata[levc]->density, 1);
     }
 
-    AMREX_D_TERM(const Real l_dx = geom[lev].CellSize(0);,
-                 const Real l_dy = geom[lev].CellSize(1);,
-                 const Real l_dz = geom[lev].CellSize(2););
+    AMREX_D_TERM(const Real l_dx = geom[levc].CellSize(0);,
+                 const Real l_dy = geom[levc].CellSize(1);,
+                 const Real l_dz = geom[levc].CellSize(2););
 
 #ifdef _OPENMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
-    for (MFIter mfi(m_leveldata[lev]->density,TilingIfNotGPU()); mfi.isValid(); ++mfi)
+    for (MFIter mfi(m_leveldata[levc]->density,TilingIfNotGPU()); mfi.isValid(); ++mfi)
     {
         Box const& bx = mfi.tilebox();
         auto const& tag = tags.array(mfi);
 
         if (tag_rho || tag_gradrho)
         {
-            Array4<Real const> const& rho = m_leveldata[lev]->density.const_array(mfi);
-            Real rhoerr = tag_rho ? rhoerr_v[lev]: std::numeric_limits<Real>::max();
-            Real gradrhoerr = tag_gradrho ? gradrhoerr_v[lev] : std::numeric_limits<Real>::max();
+            Array4<Real const> const& rho = m_leveldata[levc]->density.const_array(mfi);
+            Real rhoerr = tag_rho ? rhoerr_v[levc]: std::numeric_limits<Real>::max();
+            Real gradrhoerr = tag_gradrho ? gradrhoerr_v[levc] : std::numeric_limits<Real>::max();
             amrex::ParallelFor(bx,
             [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
             {
@@ -101,7 +101,7 @@ void incflo::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/)
             Real ylo = tag_region_lo[1];
             Real xhi = tag_region_hi[0];
             Real yhi = tag_region_hi[1];
-            auto const& problo = geom[lev].ProbLoArray();
+            auto const& problo = geom[levc].ProbLoArray();
 
 #if (AMREX_SPACEDIM == 2)
 
@@ -137,14 +137,65 @@ void incflo::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/)
             });
 #endif
         }
-    }
+    } // mfi
 
 #ifdef AMREX_USE_EB
     m_refine_cutcells = true;
     // Refine on cut cells
     if (m_refine_cutcells)
     {
-        amrex::TagCutCells(tags, m_leveldata[lev]->velocity);
+        amrex::TagCutCells(tags, m_leveldata[levc]->velocity);
     }
 #endif
+
+#ifdef INCFLO_USE_PARTICLES
+    if (m_refine_particles)
+    {
+        //
+        // This allows dynamic refinement based on the number of particles per cell
+        //
+        // Note that we must count all the particles in levels both at and above the current,
+        //      since otherwise, e.g., if the particles are all at level 1, counting particles at
+        //      level 0 will not trigger refinement when regridding so level 1 will disappear,
+        //      then come back at the next regridding
+        //
+        const auto& particles_namelist( particleData.getNames() );
+        std::unique_ptr<MultiFab> mf = std::make_unique<MultiFab>(grids[levc], dmap[levc], 1, 0);
+        mf->setVal(0.0);
+        IntVect rr = IntVect::TheUnitVector();
+        for (int lev = levc; lev <= finest_level; lev++)
+        {
+            MultiFab temp_dat(grids[lev], dmap[lev], 1, 0); temp_dat.setVal(0);
+            particleData[particles_namelist[0]]->IncrementWithTotal(temp_dat, lev);
+
+            MultiFab temp_dat_crse(grids[levc], dmap[levc], 1, 0); temp_dat_crse.setVal(0);
+
+            if (lev == levc) {
+                MultiFab::Copy(*mf, temp_dat, 0, 0, 1, 0);
+            } else {
+                for (int d = 0; d < AMREX_SPACEDIM; d++) {
+                    rr[d] *= ref_ratio[levc][d];
+                }
+                average_down(temp_dat, temp_dat_crse, 0, 1, rr);
+                MultiFab::Add(*mf, temp_dat_crse, 0, 0, 1, 0);
+            }
+        } // lev
+#ifdef _OPENMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+        for (MFIter mfi(m_leveldata[levc]->density,TilingIfNotGPU()); mfi.isValid(); ++mfi)
+        {
+            Box const& bx = mfi.tilebox();
+            auto const&  mf_arr = mf->const_array(mfi);
+            auto const& tag_arr = tags.array(mfi);
+
+            amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
+            {
+                if (mf_arr(i,j,k) > 0) {
+                    tag_arr(i,j,k) = tagval;
+                }
+            });
+        } // mfi
+    } // if m_refine_particles
+#endif
 }

src/particles/Make.package

@@ -0,0 +1,7 @@
+CEXE_sources += incflo_Tracers.cpp
+CEXE_sources += incflo_PCInit.cpp
+CEXE_sources += incflo_PCEvolve.cpp
+CEXE_sources += incflo_PCUtils.cpp
+
+CEXE_headers += ParticleData.H
+CEXE_headers += incflo_PC.H