Generalize particles (#1305)

* generalizing particle stuff * generalize particle functionality and remove unused variables * wrap with if * wrap with if
erf-model · Nov 17, 2023 · 3b5d301 · 3b5d301
1 parent 03c1754
commit 3b5d301
Show file tree

Hide file tree

Showing 46 changed files with 581 additions and 129 deletions.
diff --git a/CMake/BuildERFExe.cmake b/CMake/BuildERFExe.cmake
@@ -38,7 +38,8 @@ function(build_erf_lib erf_lib_name)
 
   if(ERF_ENABLE_PARTICLES)
     target_sources(${erf_lib_name} PRIVATE
-                   ${SRC_DIR}/Particles/TerrainFittedPC.cpp)
+                   ${SRC_DIR}/Particles/TracerPC.cpp
+                   ${SRC_DIR}/Particles/HydroPC.cpp)
     target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Particles)
     target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_PARTICLES)
   endif()

diff --git a/Docs/sphinx_doc/Particles.rst b/Docs/sphinx_doc/Particles.rst
@@ -7,9 +7,12 @@
 Particles
 =========
 
-ERF has the option to include Lagrangian particles in addition to the mesh-based solution.  Currently the
-particle functionality is very simple -- the particles are initialized one per mesh cell
-in a particular plane, and are advected by the velocity field.
+ERF has the option to include Lagrangian particles in addition to the mesh-based solution.  Currently
+there are two particle types available in ERF: tracer_particles and hydro_particles.
+The particle functionality is very simple and meant for demonstration.
+The particles are initialized one per mesh cell in a
+vertical plane at :math:`i = 3` for tracer particles and a horizontal plane at :math:`k = 23` for hydro particles.
+The tracer particles are advected by the velocity field; the hydro particles fall with a velocity determined by gravity minus drag.
 
 However, the AMReX particle data structure is very general and particles may take on a number of
 different roles in future.
@@ -32,19 +35,25 @@ One must also set
 
 ::
 
-   erf.use_tracer_particles = true
+   erf.use_tracer_particles = 1
 
-in the inputs file or on the command line at runtime.
+and / or
+
+::
 
-Currently, by default, the particles are initialized at cell centers, one per cell when the cell index is
-(3,j,k), with zero initial velocity.  They are advanced in time every time step.
+   erf.use_hydro_particles = 1
+
+in the inputs file or on the command line at runtime.
 
 Caveat: the particle information is currently output when using the AMReX-native plotfile format, but not
 when using netcdf.  Writing particles into the netcdf files is a WIP.
 
 To see an example of using the particle functionality, build the executable using gmake in Exec/DevTests/ParticlesOverWoA.
 
-To visualize the number of particles per cell as a mesh-based variable, add ``particle_count`` to the line in the inputs file
+To visualize the number of particles per cell as a mesh-based variable, add
+``tracer_particle_count`` (if you have set ``erf.use_tracer_particles``) and
+ ``hydro_particle_count`` (if you have set ``erf.use_tracer_particles``)
+to the line in the inputs file that begins
 
 ::
 

diff --git a/Exec/DevTests/ParticlesOverWoA/inputs b/Exec/DevTests/ParticlesOverWoA/inputs
@@ -24,7 +24,7 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # PARTICLES
-erf.use_tracer_particles = true
+erf.use_tracer_particles = 1
 
 # TIME STEP CONTROL
 erf.no_substepping     = 1
@@ -45,7 +45,7 @@ erf.check_int       = 10  # number of timesteps between checkpoints
 # PLOTFILES
 erf.plot_file_1     = plt     # prefix of plotfile name
 erf.plot_int_1      = 1       # number of timesteps between plotfiles
-erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens z_phys detJ dpdx dpdy pres_hse_x pres_hse_y particle_count
+erf.plot_vars_1     = density x_velocity y_velocity z_velocity pressure theta pres_hse dens_hse pert_pres pert_dens z_phys detJ dpdx dpdy pres_hse_x pres_hse_y tracer_particle_count
 
 # SOLVER CHOICE
 erf.use_gravity = true

diff --git a/Exec/RegTests/CouetteFlow/inputs_couette_x b/Exec/RegTests/CouetteFlow/inputs_couette_x
@@ -18,7 +18,6 @@ zlo.velocity    = 0.0 0.0 0.0  # for Dirichlet BC
 zhi.velocity    = 2.0 0.0 0.0  # for Dirichlet BC
 
 # TIME STEP CONTROL
-erf.use_lowM_dt     = 1
 erf.cfl             = 0.9
 
 # DIAGNOSTICS & VERBOSITY

diff --git a/Exec/RegTests/CouetteFlow/inputs_couette_y b/Exec/RegTests/CouetteFlow/inputs_couette_y
@@ -18,7 +18,6 @@ zlo.velocity    = 0.0 0.0 0.0  # for Dirichlet BC
 zhi.velocity    = 0.0 2.0 0.0  # for Dirichlet BC
 
 # TIME STEP CONTROL
-erf.use_lowM_dt     = 1
 erf.cfl             = 0.9
 
 # DIAGNOSTICS & VERBOSITY

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_WENO b/Exec/RegTests/ScalarAdvDiff/inputs_WENO
@@ -15,7 +15,6 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
-erf.use_lowM_dt    = 1
 erf.cfl            = 0.5     # cfl number for hyperbolic system
 
 erf.dycore_horiz_adv_type    = "Centered_2nd"

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z b/Exec/RegTests/ScalarAdvDiff/inputs_WENO_Z
@@ -15,7 +15,6 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
-erf.use_lowM_dt    = 1
 erf.cfl            = 0.5     # cfl number for hyperbolic system
 
 # TEST WENO SCHEMES

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_advect_shearU b/Exec/RegTests/ScalarAdvDiff/inputs_advect_shearU
@@ -15,7 +15,6 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
-erf.use_lowM_dt    = 1
 erf.cfl             = 0.9     # cfl number for hyperbolic system
 
 # DIAGNOSTICS & VERBOSITY

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_advect_shearU_msf b/Exec/RegTests/ScalarAdvDiff/inputs_advect_shearU_msf
@@ -19,7 +19,6 @@ zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
 erf.no_substepping = 1
-#erf.use_lowM_dt    = 1
 #erf.cfl             = 0.9     # cfl number for hyperbolic system
 erf.fixed_dt = 0.01
 

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_advect_shearU_no_msf b/Exec/RegTests/ScalarAdvDiff/inputs_advect_shearU_no_msf
@@ -18,7 +18,6 @@ zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
 erf.no_substepping = 1
-#erf.use_lowM_dt    = 1
 #erf.cfl             = 0.9     # cfl number for hyperbolic system
 erf.fixed_dt = 0.01
 

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_advect_uniformU b/Exec/RegTests/ScalarAdvDiff/inputs_advect_uniformU
@@ -15,7 +15,6 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
-erf.use_lowM_dt    = 1
 erf.cfl            = 0.9     # cfl number for hyperbolic system
 
 # DIAGNOSTICS & VERBOSITY

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_advect_uniformU_msf b/Exec/RegTests/ScalarAdvDiff/inputs_advect_uniformU_msf
@@ -19,7 +19,6 @@ zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
 erf.no_substepping = 1
-#erf.use_lowM_dt    = 1
 #erf.cfl            = 0.9     # cfl number for hyperbolic system
 erf.fixed_dt = 0.0005
 

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_advect_uniformU_no_msf b/Exec/RegTests/ScalarAdvDiff/inputs_advect_uniformU_no_msf
@@ -18,7 +18,6 @@ zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
 erf.no_substepping = 1
-#erf.use_lowM_dt    = 1
 #erf.cfl            = 0.9     # cfl number for hyperbolic system
 erf.fixed_dt = 0.0005
 

diff --git a/Exec/RegTests/ScalarAdvDiff/inputs_test_rayleigh b/Exec/RegTests/ScalarAdvDiff/inputs_test_rayleigh
@@ -15,7 +15,6 @@ zlo.type = "SlipWall"
 zhi.type = "SlipWall"
 
 # TIME STEP CONTROL
-erf.use_lowM_dt    = 1
 erf.cfl            = 0.9     # cfl number for hyperbolic system
 
 # DIAGNOSTICS & VERBOSITY

diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp
@@ -11,7 +11,6 @@ void
 ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
                               const amrex::MultiFab& yvel,
                               const amrex::MultiFab& cons_in,
-                              const amrex::MultiFab& cons_old,
                               amrex::MultiFab& eddyViscosity,
                               const amrex::Geometry& geom,
                               const TurbChoice& turbChoice,
@@ -373,7 +372,6 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
                                 const amrex::MultiFab& Tau11, const amrex::MultiFab& Tau22, const amrex::MultiFab& Tau33,
                                 const amrex::MultiFab& Tau12, const amrex::MultiFab& Tau13, const amrex::MultiFab& Tau23,
                                 const amrex::MultiFab& cons_in,
-                                const amrex::MultiFab& cons_old,
                                 amrex::MultiFab& eddyViscosity,
                                 amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss,
                                 const amrex::Geometry& geom,
@@ -416,7 +414,7 @@ void ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::Multi
 
     if (turbChoice.pbl_type != PBLType::None) {
         // NOTE: state_new is passed in for Cons_old (due to ptr swap in advance)
-        ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, cons_old, eddyViscosity,
+        ComputeTurbulentViscosityPBL(xvel, yvel, cons_in, eddyViscosity,
                                      geom, turbChoice, most, bc_ptr, vert_only);
     }
 }
diff --git a/Source/Diffusion/EddyViscosity.H b/Source/Diffusion/EddyViscosity.H
@@ -12,7 +12,6 @@ ComputeTurbulentViscosity (const amrex::MultiFab& xvel , const amrex::MultiFab&
                            const amrex::MultiFab& Tau11, const amrex::MultiFab& Tau22, const amrex::MultiFab& Tau33,
                            const amrex::MultiFab& Tau12, const amrex::MultiFab& Tau13, const amrex::MultiFab& Tau23,
                            const amrex::MultiFab& cons_in,
-                           const amrex::MultiFab& cons_old,
                            amrex::MultiFab& eddyViscosity,
                            amrex::MultiFab& Hfx1, amrex::MultiFab& Hfx2, amrex::MultiFab& Hfx3, amrex::MultiFab& Diss,
                            const amrex::Geometry& geom,

diff --git a/Source/Diffusion/PBLModels.cpp b/Source/Diffusion/PBLModels.cpp
@@ -19,7 +19,6 @@ void
 ComputeTurbulentViscosityPBL (const amrex::MultiFab& xvel,
                               const amrex::MultiFab& yvel,
                               const amrex::MultiFab& cons_in,
-                              const amrex::MultiFab& cons_old,
                               amrex::MultiFab& eddyViscosity,
                               const amrex::Geometry& geom,
                               const TurbChoice& turbChoice,

diff --git a/Source/ERF.H b/Source/ERF.H
@@ -27,7 +27,7 @@
 #endif
 
 #ifdef ERF_USE_PARTICLES
-#include "TerrainFittedPC.H"
+#include "ParticleData.H"
 #endif
 
 #include "prob_common.H"
@@ -66,10 +66,6 @@
 class MultiBlockContainer;
 #endif
 
-#ifdef ERF_USE_PARTICLES
-typedef amrex::ParticleContainer<AMREX_SPACEDIM, 0, 0, 0> TracerPC;
-#endif
-
 /**
  * Enum of possible coarse/fine interpolation options
 */
@@ -320,11 +316,6 @@ public:
                            int ncomp_cons = NVAR);
 #endif
 
-#ifdef ERF_USE_PARTICLES
-    std::unique_ptr<TerrainFittedPC> tracer_particles;
-    static amrex::Vector<std::string> tracer_particle_varnames;
-#endif
-
 #ifdef ERF_USE_NETCDF
     void init_from_wrfinput (int lev);
     void init_from_metgrid (int lev);
@@ -692,17 +683,19 @@ private:
                                                     ,"xvel_err", "yvel_err", "zvel_err", "pp_err"
 #endif
 #ifdef ERF_USE_PARTICLES
-                                                    ,"particle_count"
+                                                    ,"tracer_particle_count", "hydro_particle_count"
 #endif
                                                    };
 
     // algorithm choices
     static SolverChoice solverChoice;
 
-    static int verbose;
+#ifdef ERF_USE_PARTICLES
+    // Particle info
+    static ParticleData particleData;
+#endif
 
-    // flag to turn tracer particle generation on/off
-    static bool use_tracer_particles;
+    static int verbose;
 
     // Diagnostic output interval
     static int sum_interval;

diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -27,6 +27,9 @@ amrex::Real ERF::previousCPUTimeUsed = 0.0;
 Vector<AMRErrorTag> ERF::ref_tags;
 
 SolverChoice ERF::solverChoice;
+#ifdef ERF_USE_MULTIBLOCK
+ParticleData ERF::particleData;
+#endif
 
 // Time step control
 amrex::Real ERF::cfl           =  0.8;
@@ -36,12 +39,6 @@ amrex::Real ERF::init_shrink   =  1.0;
 amrex::Real ERF::change_max    =  1.1;
 int         ERF::fixed_mri_dt_ratio = 0;
 
-
-#ifdef ERF_USE_PARTICLES
-bool ERF::use_tracer_particles = false;
-amrex::Vector<std::string> ERF::tracer_particle_varnames = {AMREX_D_DECL("xvel", "yvel", "zvel")};
-#endif
-
 // Dictate verbosity in screen output
 int         ERF::verbose       = 0;
 
@@ -540,14 +537,7 @@ ERF::InitData ()
         }
 
 #ifdef ERF_USE_PARTICLES
-        // Initialize tracer particles if required
-        if (use_tracer_particles) {
-            tracer_particles = std::make_unique<TerrainFittedPC>(Geom(0), dmap[0], grids[0]);
-
-            tracer_particles->InitParticles(*z_phys_nd[0]);
-
-            Print() << "Initialized " << tracer_particles->TotalNumberOfParticles() << " tracer particles." << std::endl;
-        }
+        particleData.init_particles((amrex::ParGDBBase*)GetParGDB(),z_phys_nd);
 #endif
 
     } else { // Restart from a checkpoint
@@ -963,8 +953,8 @@ ERF::ReadParameters ()
         pp.query("fixed_mri_dt_ratio", fixed_mri_dt_ratio);
 
 #ifdef ERF_USE_PARTICLES
-        // Tracer particle toggle
-        pp.query("use_tracer_particles", use_tracer_particles);
+        particleData.init_particles((amrex::ParGDBBase*)GetParGDB(),z_phys_nd);
+        particleData.init_particle_params();
 #endif
 
 #ifdef ERF_USE_MOISTURE
@@ -1126,6 +1116,10 @@ ERF::ReadParameters ()
     solverChoice.pp_prefix = pp_prefix;
 #endif
 
+#ifdef ERF_USE_MULTIBLOCK
+    Particleinit_particle_params();
+#endif
+
     solverChoice.init_params(max_level);
     if (verbose > 0) {
         solverChoice.display();

diff --git a/Source/IO/Checkpoint.cpp b/Source/IO/Checkpoint.cpp
@@ -173,9 +173,7 @@ ERF::WriteCheckpointFile () const
    }
 
 #ifdef ERF_USE_PARTICLES
-   if (use_tracer_particles) {
-       tracer_particles->Checkpoint(checkpointname, "tracers", true, tracer_particle_varnames);
-   }
+   particleData.Checkpoint(checkpointname);
 #endif
 
 #ifdef ERF_USE_NETCDF
@@ -378,11 +376,7 @@ ERF::ReadCheckpointFile ()
     }
 
 #ifdef ERF_USE_PARTICLES
-   if (use_tracer_particles) {
-       tracer_particles = std::make_unique<TerrainFittedPC>(Geom(0), dmap[0], grids[0]);
-       std::string tracer_file("tracers");
-       tracer_particles->Restart(restart_chkfile, tracer_file);
-   }
+   particleData.Restart((amrex::ParGDBBase*)GetParGDB(),restart_chkfile);
 #endif
 
 #ifdef ERF_USE_NETCDF

diff --git a/Source/IO/ERF_WriteScalarProfiles.cpp b/Source/IO/ERF_WriteScalarProfiles.cpp
@@ -145,7 +145,7 @@ ERF::sum_integrated_quantities(Real time)
  * @param mf MultiFab from which we wish to sample data
  */
 void
-ERF::sample_points(int lev, Real time, IntVect cell, MultiFab& mf)
+ERF::sample_points(int /*lev*/, Real time, IntVect cell, MultiFab& mf)
 {
     int datwidth = 14;
 

diff --git a/Source/IO/Plotfile.cpp b/Source/IO/Plotfile.cpp
@@ -613,11 +613,19 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
 #endif
 
 #ifdef ERF_USE_PARTICLES
-        if (containerHasElement(plot_var_names, "particle_count"))
+        if (containerHasElement(plot_var_names, "tracer_particle_count"))
         {
             MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
             temp_dat.setVal(0);
-            tracer_particles->Increment(temp_dat, lev);
+            particleData.tracer_particles->Increment(temp_dat, lev);
+            MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
+            mf_comp += 1;
+        }
+        if (containerHasElement(plot_var_names, "hydro_particle_count"))
+        {
+            MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
+            temp_dat.setVal(0);
+            particleData.hydro_particles->Increment(temp_dat, lev);
             MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
             mf_comp += 1;
         }
@@ -827,9 +835,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             writeJobInfo(plotfilename);
 
 #ifdef ERF_USE_PARTICLES
-            if (use_tracer_particles) {
-                tracer_particles->Checkpoint(plotfilename, "tracers", true, tracer_particle_varnames);
-            }
+            particleData.Checkpoint(plotfilename);
 #endif
 #ifdef ERF_USE_HDF5
         } else if (plotfile_type == "hdf5" || plotfile_type == "HDF5") {
@@ -917,9 +923,7 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
             writeJobInfo(plotfilename);
 
 #ifdef ERF_USE_PARTICLES
-            if (use_tracer_particles) {
-                tracer_particles->Checkpoint(plotfilename, "tracers", true, tracer_particle_varnames);
-            }
+            particleData.Checkpoint(plotfilename);
 #endif
 #ifdef ERF_USE_NETCDF
         } else if (plotfile_type == "netcdf" || plotfile_type == "NetCDF") {