Initial implementation of a generalized ERF particle interface

Source/ERF.H

@@ -26,10 +26,6 @@
 #include <AMReX_MemProfiler.H>
 #endif
 
-#ifdef ERF_USE_PARTICLES
-#include "ParticleData.H"
-#endif
-
 #include "prob_common.H"
 
 #include <IndexDefines.H>
@@ -44,6 +40,10 @@
 #include <ERF_PhysBCFunct.H>
 #include <ERF_FillPatcher.H>
 
+#ifdef ERF_USE_PARTICLES
+#include "ParticleData.H"
+#endif
+
 #include "Microphysics.H"
 
 #ifdef ERF_USE_RRTMGP
@@ -660,18 +660,32 @@ private:
  ,"qt", "qp", "qv", "qc", "qi", "qrain", "qsnow", "qgraup"
 #ifdef ERF_COMPUTE_ERROR
  ,"xvel_err", "yvel_err", "zvel_err", "pp_err"
-#endif
-#ifdef ERF_USE_PARTICLES
- ,"tracer_particle_count", "hydro_particle_count"
 #endif
  };
 
  // algorithm choices
  static SolverChoice solverChoice;
 
 #ifdef ERF_USE_PARTICLES
- // Particle info
+ // Particle container with all particle species
  static ParticleData particleData;
+
+ // variables and functions for tracers particles
+ bool m_use_tracer_particles; /*!< tracer particles that advect with flow */
+ bool m_use_hydro_particles; /*!< tracer particles that fall with gravity */
+
+ /*! Read tracer and hydro particles parameters */
+ void readTracersParams();
+
+ /*! Initialize tracer and hydro particles */
+ void initializeTracers ( amrex::ParGDBBase*,
+ const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& );
+ /*! Evolve tracers and hydro particles */
+ void evolveTracers( int,
+ amrex::Real,
+ amrex::Vector<amrex::Vector<amrex::MultiFab>>&,
+ const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& );
+
 #endif
 
  static int verbose;

Source/ERF.cpp

@@ -537,7 +537,7 @@ ERF::InitData ()
  }
 
 #ifdef ERF_USE_PARTICLES
- particleData.init_particles((amrex::ParGDBBase*)GetParGDB(),z_phys_nd);
+ initializeTracers((amrex::ParGDBBase*)GetParGDB(),z_phys_nd);
 #endif
 
  } else { // Restart from a checkpoint
@@ -945,10 +945,6 @@ ERF::ReadParameters ()
  pp.query("fixed_fast_dt", fixed_fast_dt);
  pp.query("fixed_mri_dt_ratio", fixed_mri_dt_ratio);
 
-#ifdef ERF_USE_PARTICLES
- particleData.init_particle_params();
-#endif
-
  // If this is set, it must be even
  if (fixed_mri_dt_ratio > 0 && (fixed_mri_dt_ratio%2 != 0) )
  {
@@ -1094,12 +1090,12 @@ ERF::ReadParameters ()
  if (!iterate) SetIterateToFalse();
  }
 
-#ifdef ERF_USE_MULTIBLOCK
- solverChoice.pp_prefix = pp_prefix;
+#ifdef ERF_USE_PARTICLES
+ readTracersParams();
 #endif
 
-#ifdef ERF_USE_PARTICLES
- particleData.init_particle_params();
+#ifdef ERF_USE_MULTIBLOCK
+ solverChoice.pp_prefix = pp_prefix;
 #endif
 
  solverChoice.init_params(max_level);

Source/IO/Plotfile.cpp

@@ -63,17 +63,20 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
  }
  for (int i = 0; i < derived_names.size(); ++i) {
  if ( containerHasElement(plot_var_names, derived_names[i]) ) {
-#ifdef ERF_USE_PARTICLES
- if ( (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) &&
- (particleData.use_tracer_particles || (derived_names[i] != "tracer_particle_count") ) &&
- (particleData.use_hydro_particles || (derived_names[i] != "hydro_particle_count") ) ) {
-#else
  if (solverChoice.use_terrain || (derived_names[i] != "z_phys" && derived_names[i] != "detJ") ) {
-#endif
  tmp_plot_names.push_back(derived_names[i]);
  }
  }
  }
+#ifdef ERF_USE_PARTICLES
+ const auto& particles_namelist( particleData.getNamesUnalloc() );
+ for (auto it = particles_namelist.cbegin(); it != particles_namelist.cend(); ++it) {
+ std::string tmp( (*it)+"_count" );
+ if (containerHasElement(plot_var_names, tmp) ) {
+ tmp_plot_names.push_back(tmp);
+ }
+ }
+#endif
 
  // Check to see if we found all the requested variables
  for (const auto& plot_name : plot_var_names) {
@@ -691,21 +694,15 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
  }
 
 #ifdef ERF_USE_PARTICLES
- if (containerHasElement(plot_var_names, "tracer_particle_count"))
- {
- MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
- temp_dat.setVal(0);
- 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;
+ const auto& particles_namelist( particleData.getNames() );
+ for (ParticlesNamesVector::size_type i = 0; i < particles_namelist.size(); i++) {
+ if (containerHasElement(plot_var_names, std::string(particles_namelist[i]+"_count"))) {
+ MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 0);
+ temp_dat.setVal(0);
+ particleData[particles_namelist[i]]->Increment(temp_dat, lev);
+ MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
+ mf_comp += 1;
+ }
  }
 #endif
 

Source/Particles/ERFPC.H

@@ -0,0 +1,173 @@
+#ifndef ERF_PC_H_
+#define ERF_PC_H_
+
+#ifdef ERF_USE_PARTICLES
+
+#include <string>
+#include <AMReX_Particles.H>
+
+struct ERFParticlesRealIdx
+{
+ enum {
+ vx = 0,
+ vy,
+ vz,
+ mass,
+ ncomps
+ };
+};
+
+struct ERFParticlesIntIdx
+{
+ enum {
+ i = 0,
+ j,
+ k,
+ ncomps
+ };
+};
+
+namespace ERFParticleInitializations
+{
+ /* list of particle initializations */
+ const std::string init_default = "default";
+ const std::string init_uniform = "uniform";
+
+}
+
+namespace ERFParticleNames
+{
+ const std::string tracers = "tracer_particles";
+ const std::string hydro = "hydro_particles";
+}
+
+struct ERFParticlesAssignor
+{
+ template <typename P>
+ AMREX_GPU_HOST_DEVICE
+ amrex::IntVect operator() ( P const& p,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
+ const amrex::Box& domain ) const noexcept
+ {
+ /* TODO: What is this about? */
+ amrex::IntVect iv(
+ AMREX_D_DECL( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
+ int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
+ p.idata(ERFParticlesIntIdx::k) ) );
+ iv[0] += domain.smallEnd()[0];
+ iv[1] += domain.smallEnd()[1];
+ return iv;
+ }
+};
+
+class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdx::ncomps,
+ ERFParticlesIntIdx::ncomps,
+ 0,
+ 0,
+ amrex::DefaultAllocator,
+ ERFParticlesAssignor >
+{
+ public:
+
+ /*! Constructor */
+ ERFPC ( amrex::ParGDBBase* a_gdb,
+ const std::string& a_name = "particles" )
+ : amrex::ParticleContainer< ERFParticlesRealIdx::ncomps,
+ ERFParticlesIntIdx::ncomps,
+ 0,
+ 0,
+ amrex::DefaultAllocator,
+ ERFParticlesAssignor> (a_gdb)
+ {
+ BL_PROFILE("ERFPCPC::ERFPC()");
+ m_name = a_name;
+ readInputs();
+ }
+
+ /*! Constructor */
+ ERFPC ( const amrex::Geometry& a_geom,
+ const amrex::DistributionMapping& a_dmap,
+ const amrex::BoxArray& a_ba,
+ const std::string& a_name = "particles" )
+ : amrex::ParticleContainer< ERFParticlesRealIdx::ncomps,
+ ERFParticlesIntIdx::ncomps,
+ 0,
+ 0,
+ amrex::DefaultAllocator,
+ ERFParticlesAssignor> ( a_geom, a_dmap, a_ba )
+ {
+ BL_PROFILE("ERFPCPC::ERFPC()");
+ m_name = a_name;
+ readInputs();
+ }
+
+ /*! Initialize particles in domain */
+ virtual void InitializeParticles(const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);
+
+ /*! Evolve particles for one time step */
+ virtual void EvolveParticles( int,
+ amrex::Real,
+ amrex::Vector<amrex::Vector<amrex::MultiFab>>&,
+ const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& );
+
+ /*! Get real-type particle attribute names */
+ virtual amrex::Vector<std::string> varNames() const
+ {
+ BL_PROFILE("ERFPCPC::varNames()");
+ return {AMREX_D_DECL("xvel","yvel","zvel"),"mass"};
+ }
+
+ /*! Specify if particles should advect with flow */
+ inline void setAdvectWithFlow (bool a_flag)
+ {
+ BL_PROFILE("ERFPCPC::setAdvectWithFlow()");
+ m_advect_w_flow = a_flag;
+ }
+ /*! Specify if particles fall under gravity */
+ inline void setAdvectWithGravity (bool a_flag)
+ {
+ BL_PROFILE("ERFPCPC::setAdvectWithGravity()");
+ m_advect_w_gravity = a_flag;
+ }
+
+ protected:
+
+ bool m_advect_w_flow; /*!< advect with flow velocity */
+ bool m_advect_w_gravity; /*!< advect under gravitational force */
+
+ std::string m_name; /*!< name of this particle species */
+
+ std::string m_initialization_type; /*!< initial particle distribution type */
+ int m_ppc_init; /*!< initial number of particles per cell */
+
+ /*! read inputs from file */
+ virtual void readInputs();
+
+ /*! Default particle initialization */
+ void initializeParticlesUniformDistribution(const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);
+
+ /*! Uses midpoint method to advance particles using flow velocity. */
+ void AdvectWithFlow( amrex::MultiFab*,
+ int,
+ amrex::Real,
+ const std::unique_ptr<amrex::MultiFab>& );
+
+ /*! Uses midpoint method to advance particles falling under gravity. */
+ void AdvectWithGravity( int,
+ amrex::Real,
+ const std::unique_ptr<amrex::MultiFab>& );
+
+ private:
+
+ /*! Default particle initialization */
+ void initializeParticlesDefault(const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);
+ /*! Default initialization for tracer particles for WoA case (ref: AA) */
+ void initializeParticlesDefaultTracersWoA(const std::unique_ptr<amrex::MultiFab>& a_ptr=nullptr);
+ /*! Default initialization for hydro particles (ref: AA) */
+ void initializeParticlesDefaultHydro(const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);
+
+};
+
+#endif
+#endif