Skip to content

Commit

Permalink
initial implementation of a generalized ERF particle interface; retai…
Browse files Browse the repository at this point in the history
…ned tracer and hydro particles as specific types; added uniform distribution of particles as an initialization option
  • Loading branch information
debog committed Jan 4, 2024
1 parent 53179f4 commit 9d432d6
Show file tree
Hide file tree
Showing 14 changed files with 1,117 additions and 753 deletions.
30 changes: 22 additions & 8 deletions Source/ERF.H
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,6 @@
#include <AMReX_MemProfiler.H>
#endif

#ifdef ERF_USE_PARTICLES
#include "ParticleData.H"
#endif

#include "prob_common.H"

#include <IndexDefines.H>
Expand All @@ -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
Expand Down Expand Up @@ -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;
Expand Down
14 changes: 5 additions & 9 deletions Source/ERF.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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) )
{
Expand Down Expand Up @@ -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);
Expand Down
39 changes: 18 additions & 21 deletions Source/IO/Plotfile.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down Expand Up @@ -667,21 +670,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

Expand Down
173 changes: 173 additions & 0 deletions Source/Particles/ERFPC.H
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 9d432d6

Please sign in to comment.