Generalized particle restart to allow restarting other particle conta…

…iners (#1862) * generalized the particle restart * separated defining the particle container and initializing particles to avoid initializing new particles for a restart simulation * created an interface to restart Lagrangian microphysics particle containers * added missing header
erf-model · Oct 8, 2024 · 13a7ce3 · 13a7ce3
1 parent 1d497e7
commit 13a7ce3
Show file tree

Hide file tree

Showing 7 changed files with 72 additions and 15 deletions.
diff --git a/Source/ERF.H b/Source/ERF.H
@@ -930,6 +930,10 @@ private:
     /*! Initialize tracer and hydro particles */
     void initializeTracers ( amrex::ParGDBBase*,
                              const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& );
+
+    /*! Restart tracer and hydro particles */
+    void restartTracers ( amrex::ParGDBBase*, const std::string& );
+
     /*! Evolve tracers and hydro particles */
     void evolveTracers( int,
                         amrex::Real,

diff --git a/Source/ERF.cpp b/Source/ERF.cpp
@@ -722,6 +722,14 @@ ERF::InitData_post ()
             }
         }
 
+#ifdef ERF_USE_PARTICLES
+        if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
+            for (int lev = 0; lev <= finest_level; lev++) {
+                dynamic_cast<LagrangianMicrophysics&>(*micro).initParticles(z_phys_nd[lev]);
+            }
+        }
+#endif
+
     } else { // Restart from a checkpoint
 
         restart();

diff --git a/Source/IO/ERF_Checkpoint.cpp b/Source/IO/ERF_Checkpoint.cpp
@@ -454,7 +454,10 @@ ERF::ReadCheckpointFile ()
     }
 
 #ifdef ERF_USE_PARTICLES
-   particleData.Restart((ParGDBBase*)GetParGDB(),restart_chkfile);
+    restartTracers((ParGDBBase*)GetParGDB(),restart_chkfile);
+    if (Microphysics::modelType(solverChoice.moisture_type) == MoistureModelType::Lagrangian) {
+        dynamic_cast<LagrangianMicrophysics&>(*micro).restartParticles((ParGDBBase*)GetParGDB(),restart_chkfile);
+    }
 #endif
 
 #ifdef ERF_USE_NETCDF

diff --git a/Source/Microphysics/ERF_LagrangianMicrophysics.H b/Source/Microphysics/ERF_LagrangianMicrophysics.H
@@ -111,6 +111,18 @@ public:
         return m_moist_model->Qstate_Size();
     }
 
+    /*! \brief initialize particles in particle container */
+    inline void initParticles ( std::unique_ptr<amrex::MultiFab>& z_phys_nd /*!< Nodal z heights */ )
+    {
+        m_moist_model->InitParticles(z_phys_nd);
+    }
+
+    /*! \brief restart particles in particle container */
+    inline void restartParticles ( amrex::ParGDBBase* a_gdb, const std::string& a_fname)
+    {
+        m_moist_model->RestartParticles(a_gdb, a_fname);
+    }
+
     /*! \brief get the particle container from the moisture model */
     inline ERFPC* getParticleContainer () const
     {

diff --git a/Source/Microphysics/Null/ERF_NullMoistLagrangian.H b/Source/Microphysics/Null/ERF_NullMoistLagrangian.H
@@ -6,6 +6,7 @@
 
 #ifdef ERF_USE_PARTICLES
 
+#include <AMReX_AmrParGDB.H>
 #include "ERF_NullMoist.H"
 
 /* forward declaration */
@@ -24,6 +25,14 @@ public:
     /*! \brief Default destructor */
     virtual ~NullMoistLagrangian () = default;
 
+    /*! \brief Initialize particles */
+    virtual
+    void InitParticles ( std::unique_ptr<amrex::MultiFab>& ) { }
+
+    /*! \brief Restart particles */
+    virtual
+    void RestartParticles ( amrex::ParGDBBase*, const std::string& ) { }
+
     /*! \brief get the particle container */
     virtual ERFPC* getParticleContainer ()
     {

diff --git a/Source/Particles/ERFTracers.cpp b/Source/Particles/ERFTracers.cpp
@@ -62,6 +62,41 @@ void ERF::initializeTracers ( ParGDBBase* a_gdb,
     return;
 }
 
+/*! Restart tracer and hydro particles */
+void ERF::restartTracers ( ParGDBBase* a_gdb,
+                           const std::string& a_fname )
+{
+    auto& namelist_unalloc( particleData.getNamesUnalloc() );
+
+    for (auto it = namelist_unalloc.begin(); it != namelist_unalloc.end(); ++it) {
+
+        std::string species_name( *it );
+
+        if (species_name == ERFParticleNames::tracers) {
+
+            AMREX_ASSERT(m_use_tracer_particles);
+            ERFPC* pc = new ERFPC(a_gdb, ERFParticleNames::tracers);
+            pc->Restart(a_fname, ERFParticleNames::tracers);
+            amrex::Print() << "Restarted " << pc->TotalNumberOfParticles() << " tracer particles.\n";
+            particleData.pushBack(ERFParticleNames::tracers, pc);
+
+        } else if (species_name == ERFParticleNames::hydro) {
+
+            AMREX_ASSERT(m_use_hydro_particles);
+            ERFPC* pc = new ERFPC(a_gdb, ERFParticleNames::hydro);
+            pc->Restart(a_fname, ERFParticleNames::hydro);
+            amrex::Print() << "Restarted " << pc->TotalNumberOfParticles() << " hydro particles.\n";
+            particleData.pushBack(ERFParticleNames::hydro, pc);
+
+        }
+    }
+
+    if (m_use_tracer_particles) namelist_unalloc.remove( ERFParticleNames::tracers );
+    if (m_use_hydro_particles)  namelist_unalloc.remove( ERFParticleNames::hydro );
+
+    return;
+}
+
 /*! Evolve tracers and hydro particles for one time step*/
 void ERF::evolveTracers ( int                                        a_lev,
                           Real                                       a_dt_lev,

diff --git a/Source/Particles/ERF_ParticleData.H b/Source/Particles/ERF_ParticleData.H
@@ -79,20 +79,6 @@ class ParticleData
             }
         }
 
-        /*! Read from restart file */
-        void Restart ( amrex::ParGDBBase* a_gdb, const std::string& a_fname )
-        {
-            BL_PROFILE("ParticleData::Restart()");
-            AMREX_ASSERT(isEmpty());
-            for (auto it = m_namelist_unalloc.begin(); it != m_namelist_unalloc.end(); ++it) {
-                std::string species_name( *it );
-                ERFPC* pc = new ERFPC( a_gdb, species_name );
-                pc->Restart(a_fname, species_name );
-                pushBack( species_name, pc );
-            }
-            m_namelist_unalloc.clear();
-        }
-
         /*! Get mesh plot quantities from each particle container */
         void GetMeshPlotVarNames ( amrex::Vector<std::string>& a_names ) const
         {