Implemented plotting Eulerian quantities derived from particle contai…

…ners; added a mass density function to tracer particle container as an example (erf-model#1463)

CMake/BuildERFExe.cmake

@@ -38,6 +38,7 @@ function(build_erf_lib erf_lib_name)
     target_sources(${erf_lib_name} PRIVATE
                    ${SRC_DIR}/Particles/ERFPCEvolve.cpp
                    ${SRC_DIR}/Particles/ERFPCInitializations.cpp
+                   ${SRC_DIR}/Particles/ERFPCUtils.cpp
                    ${SRC_DIR}/Particles/ERFTracers.cpp)
     target_include_directories(${erf_lib_name} PUBLIC ${SRC_DIR}/Particles)
     target_compile_definitions(${erf_lib_name} PUBLIC ERF_USE_PARTICLES)

Exec/DevTests/ParticlesOverWoA/inputs

@@ -44,7 +44,7 @@ erf.check_int       = 10  # number of timesteps between checkpoints
 # PLOTFILES
 erf.plot_file_1     = plt     # prefix of plotfile name
 erf.plot_int_1      = 5       # 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 tracer_particles_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_particles_count tracer_particles_mass_density
 
 # SOLVER CHOICE
 erf.use_gravity = true

Source/ERF.H

@@ -432,6 +432,8 @@ private:
 
     // set which variables and derived quantities go into plotfiles
     void setPlotVariables (const std::string& pp_plot_var_names, amrex::Vector<std::string>& plot_var_names);
+    // append variables to plot
+    void appendPlotVariables (const std::string& pp_plot_var_names, amrex::Vector<std::string>& plot_var_names);
 
 #ifdef ERF_USE_NETCDF
     //! Write plotfile using NETCDF

Source/ERF.cpp

@@ -778,6 +778,10 @@ ERF::InitData ()
         for (int lev = 0; lev <= finest_level; ++lev) micro.Update_Micro_Vars_Lev(lev, vars_new[lev][Vars::cons]);
     }
 
+    // check for additional plotting variables that are available after particle containers
+    // are setup.
+    const std::string& pv1 = "plot_vars_1"; appendPlotVariables(pv1,plot_var_names_1);
+    const std::string& pv2 = "plot_vars_2"; appendPlotVariables(pv2,plot_var_names_2);
 
     if ( restart_chkfile.empty() && (m_check_int > 0 || m_check_per > 0.) )
     {

Source/IO/Plotfile.cpp

@@ -77,15 +77,52 @@ ERF::setPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>
     }
 #endif
 
-    // Check to see if we found all the requested variables
+    plot_var_names = tmp_plot_names;
+}
+
+void
+ERF::appendPlotVariables (const std::string& pp_plot_var_names, Vector<std::string>& a_plot_var_names)
+{
+    ParmParse pp(pp_prefix);
+
+    Vector<std::string> plot_var_names(0);
+    if (pp.contains(pp_plot_var_names.c_str())) {
+        std::string nm;
+        int nPltVars = pp.countval(pp_plot_var_names.c_str());
+        for (int i = 0; i < nPltVars; i++) {
+            pp.get(pp_plot_var_names.c_str(), nm, i);
+            // Add the named variable to our list of plot variables
+            // if it is not already in the list
+            if (!containerHasElement(plot_var_names, nm)) {
+                plot_var_names.push_back(nm);
+            }
+        }
+    }
+
+    Vector<std::string> tmp_plot_names(0);
+#ifdef ERF_USE_PARTICLES
+    Vector<std::string> particle_mesh_plot_names;
+    particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
+    for (int i = 0; i < particle_mesh_plot_names.size(); i++) {
+        std::string tmp(particle_mesh_plot_names[i]);
+        if (containerHasElement(plot_var_names, tmp) ) {
+            tmp_plot_names.push_back(tmp);
+        }
+    }
+#endif
+
+    for (int i = 0; i < tmp_plot_names.size(); i++) {
+        a_plot_var_names.push_back( tmp_plot_names[i] );
+    }
+
+    // Finally, check to see if we found all the requested variables
     for (const auto& plot_name : plot_var_names) {
-      if (!containerHasElement(tmp_plot_names, plot_name)) {
-           if (amrex::ParallelDescriptor::IOProcessor()) {
-               Warning("\nWARNING: Requested to plot variable '" + plot_name + "' but it is not available");
-           }
-      }
+        if (!containerHasElement(a_plot_var_names, plot_name)) {
+             if (amrex::ParallelDescriptor::IOProcessor()) {
+                 Warning("\nWARNING: Requested to plot variable '" + plot_name + "' but it is not available");
+             }
+        }
     }
-    plot_var_names = tmp_plot_names;
 }
 
 // set plotfile variable names
@@ -722,6 +759,19 @@ ERF::WritePlotFile (int which, Vector<std::string> plot_var_names)
                 mf_comp += 1;
             }
         }
+
+        Vector<std::string> particle_mesh_plot_names(0);
+        particleData.GetMeshPlotVarNames( particle_mesh_plot_names );
+        for (int i = 0; i < particle_mesh_plot_names.size(); i++) {
+            std::string plot_var_name(particle_mesh_plot_names[i]);
+            if (containerHasElement(plot_var_names, plot_var_name) ) {
+                MultiFab temp_dat(mf[lev].boxArray(), mf[lev].DistributionMap(), 1, 1);
+                temp_dat.setVal(0);
+                particleData.GetMeshPlotVar(plot_var_name, temp_dat, lev);
+                MultiFab::Copy(mf[lev], temp_dat, 0, mf_comp, 1, 0);
+                mf_comp += 1;
+            }
+        }
 #endif
 
 #ifdef ERF_COMPUTE_ERROR

Source/Particles/ERFPC.H

@@ -160,6 +160,13 @@ class ERFPC : public amrex::ParticleContainer<  ERFParticlesRealIdxAoS::ncomps,
             return {AMREX_D_DECL("xvel","yvel","zvel"),"mass"};
         }
 
+        /*! Get real-type particle attribute names */
+        virtual amrex::Vector<std::string> meshPlotVarNames () const
+        {
+            BL_PROFILE("ERFPCPC::varNames()");
+            return {"mass_density"};
+        }
+
         /*! Uses midpoint method to advance particles using flow velocity. */
         virtual void AdvectWithFlow (  amrex::MultiFab*,
                                        int,
@@ -171,6 +178,21 @@ class ERFPC : public amrex::ParticleContainer<  ERFParticlesRealIdxAoS::ncomps,
                                          amrex::Real,
                                          const std::unique_ptr<amrex::MultiFab>& );
 
+        /*! Compute mass density */
+        virtual void massDensity ( amrex::MultiFab&, const int&, const int& a_comp = 0) const;
+
+        /*! Compute mesh variable from particles */
+        virtual void computeMeshVar(    const std::string&  a_var_name,
+                                        amrex::MultiFab&    a_mf,
+                                        const int           a_lev) const
+        {
+            if (a_var_name == "mass_density") {
+                massDensity( a_mf, a_lev );
+            } else {
+                a_mf.setVal(0.0);
+            }
+        }
+
         /*! Specify if particles should advect with flow */
         inline void setAdvectWithFlow (bool a_flag)
         {

Source/Particles/ERFPCEvolve.cpp

@@ -24,6 +24,8 @@ void ERFPC::EvolveParticles ( int                                        a_lev,
     if (m_advect_w_gravity) {
         AdvectWithGravity( a_lev, a_dt_lev, a_z_phys_nd[a_lev] );
     }
+
+    Redistribute();
     return;
 }
 

Source/Particles/ERFPCInitializations.cpp

@@ -149,7 +149,7 @@ void ERFPC::initializeParticlesDefaultTracersWoA (const std::unique_ptr<amrex::M
                     vy_ptr[offset_arr(i,j,k)] = v[1];
                     vz_ptr[offset_arr(i,j,k)] = v[2];
 
-                    mass_ptr[offset_arr(i,j,k)] = 0.0;
+                    mass_ptr[offset_arr(i,j,k)] = 1.0e-6;
                }
             });
 
@@ -178,7 +178,7 @@ void ERFPC::initializeParticlesDefaultTracersWoA (const std::unique_ptr<amrex::M
                     vy_ptr[offset_arr(i,j,k)] = v[1];
                     vz_ptr[offset_arr(i,j,k)] = v[2];
 
-                    mass_ptr[offset_arr(i,j,k)] = 0.0;
+                    mass_ptr[offset_arr(i,j,k)] = 1.0e-6;
                }
             });
 
@@ -281,7 +281,7 @@ void ERFPC::initializeParticlesDefaultHydro (const std::unique_ptr<amrex::MultiF
                     vy_ptr[offset_arr(i,j,k)] = v[1];
                     vz_ptr[offset_arr(i,j,k)] = v[2];
 
-                    mass_ptr[offset_arr(i,j,k)] = 0.0;
+                    mass_ptr[offset_arr(i,j,k)] = 1.0e-6;
                }
             });
 
@@ -310,7 +310,7 @@ void ERFPC::initializeParticlesDefaultHydro (const std::unique_ptr<amrex::MultiF
                     vy_ptr[offset_arr(i,j,k)] = v[1];
                     vz_ptr[offset_arr(i,j,k)] = v[2];
 
-                    mass_ptr[offset_arr(i,j,k)] = 0.0;
+                    mass_ptr[offset_arr(i,j,k)] = 1.0e-6;
                }
             });
 
@@ -415,7 +415,7 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr<amrex:
                     vy_ptr[n] = v[1];
                     vz_ptr[n] = v[2];
 
-                    mass_ptr[n] = 0.0;
+                    mass_ptr[n] = 1.0e-6;
                }
             });
 
@@ -446,7 +446,7 @@ void ERFPC::initializeParticlesUniformDistribution (const std::unique_ptr<amrex:
                     vy_ptr[n] = v[1];
                     vz_ptr[n] = v[2];
 
-                    mass_ptr[n] = 0.0;
+                    mass_ptr[n] = 1.0e-6;
                }
             });
         }

Source/Particles/ERFPCUtils.cpp

@@ -0,0 +1,42 @@
+#include <AMReX_ParticleInterpolators.H>
+#include <ERF_Constants.H>
+#include <ERFPC.H>
+
+#ifdef ERF_USE_PARTICLES
+
+using namespace amrex;
+
+void ERFPC::massDensity ( MultiFab&  a_mf,
+                          const int& a_lev,
+                          const int& a_comp ) const
+{
+    BL_PROFILE("ERFPC::massDensity()");
+
+    AMREX_ASSERT(OK());
+    AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0);
+
+    const auto& geom = Geom(a_lev);
+    const auto plo = geom.ProbLoArray();
+    const auto dxi = geom.InvCellSizeArray();
+
+    const Real inv_cell_volume = dxi[0]*dxi[1]*dxi[2];
+    a_mf.setVal(0.0);
+
+    ParticleToMesh( *this, a_mf, a_lev,
+        [=] AMREX_GPU_DEVICE (  const ERFPC::ParticleTileType::ConstParticleTileDataType& ptd,
+                                int i, Array4<Real> const& rho)
+        {
+            auto p = ptd.m_aos[i];
+            ParticleInterpolator::Linear interp(p, plo, dxi);
+            interp.ParticleToMesh ( p, rho, 0, a_comp, 1,
+                [=] AMREX_GPU_DEVICE ( const ERFPC::ParticleType&, int)
+                {
+                    auto mass = ptd.m_rdata[ERFParticlesRealIdxSoA::mass][i];
+                    return mass*inv_cell_volume;
+                });
+        });
+
+    return;
+}
+
+#endif

Source/Particles/Make.package

@@ -1,6 +1,7 @@
 CEXE_sources += ERFTracers.cpp
 CEXE_sources += ERFPCInitializations.cpp
 CEXE_sources += ERFPCEvolve.cpp
+CEXE_sources += ERFPCUtils.cpp
 
 CEXE_headers += ParticleData.H
 CEXE_headers += ERFPC.H

Source/Particles/ParticleData.H

@@ -11,6 +11,7 @@
 
 #include <AMReX_ParmParse.H>
 #include <AMReX_Print.H>
+#include <AMReX_Vector.H>
 #include <AMReX_Gpu.H>
 
 #include <ERFPC.H>
@@ -74,6 +75,45 @@ class ParticleData
             m_namelist_unalloc.clear();
         }
 
+        /*! Get mesh plot quantities from each particle container */
+        void GetMeshPlotVarNames ( amrex::Vector<std::string>& a_names ) const
+        {
+            BL_PROFILE("ParticleData::GetMeshPlotVarNames()");
+            a_names.clear();
+            for (ParticlesNamesVector::size_type i = 0; i < m_namelist.size(); i++) {
+                auto name( m_namelist[i] );
+                auto particles( m_particle_species.at(name) );
+
+                auto var_names = particles->meshPlotVarNames();
+                for (int n = 0; n < var_names.size(); n++) {
+                    a_names.push_back( std::string(name+"_"+var_names[n]) );
+                }
+            }
+        }
+
+        void GetMeshPlotVar (   const std::string& a_var_name,
+                                amrex::MultiFab&   a_mf,
+                                const int          a_lev )
+        {
+            BL_PROFILE("ParticleData::GetMeshPlotVar()");
+            for (ParticlesNamesVector::size_type i = 0; i < m_namelist.size(); i++) {
+                auto particle_name( m_namelist[i] );
+                auto particles( m_particle_species.at(particle_name) );
+
+                auto particle_var_names = particles->meshPlotVarNames();
+
+                for (int n = 0; n < particle_var_names.size(); n++) {
+
+                    std::string var_name = particle_name+"_"+particle_var_names[n];
+                    if ( var_name == a_var_name ) {
+                        particles->computeMeshVar(particle_var_names[n], a_mf, a_lev);
+                        return;
+                    }
+                }
+            }
+            amrex::Abort("Requested var_name not found in ParticleData::GetMeshPlotVar");
+        }
+
         /*! Redistribute/rebalance particles data */
         inline void Redistribute ()
         {