ERF
Energy Research and Forecasting: An Atmospheric Modeling Code
ERFPC.H
Go to the documentation of this file.
1 #ifndef ERF_PC_H_
2 #define ERF_PC_H_
3 
4 #ifdef ERF_USE_PARTICLES
5 
6 #include <string>
7 #include <AMReX_Particles.H>
8 
9 struct ERFParticlesIntIdxAoS
10 {
11  enum {
12  k = 0,
13  ncomps
14  };
15 };
16 
17 struct ERFParticlesRealIdxAoS
18 {
19  enum {
20  ncomps = 0
21  };
22 };
23 
24 struct ERFParticlesIntIdxSoA
25 {
26  enum {
27  ncomps = 0
28  };
29 };
30 
31 struct ERFParticlesRealIdxSoA
32 {
33  enum {
34  vx = 0,
35  vy,
36  vz,
37  mass,
38  ncomps
39  };
40 };
41 
42 namespace ERFParticleInitializations
43 {
44  /* list of particle initializations */
45  const std::string init_box_uniform = "box";
46 }
47 
48 namespace ERFParticleNames
49 {
50  const std::string tracers = "tracer_particles";
51  const std::string hydro = "hydro_particles";
52 }
53 
54 struct ERFParticlesAssignor
55 {
56  template <typename P>
57  AMREX_GPU_HOST_DEVICE
58  amrex::IntVect operator() ( P const& p,
59  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
60  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi,
61  const amrex::Box& domain ) const noexcept
62  {
63  amrex::IntVect iv(
64  AMREX_D_DECL( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])),
65  int(amrex::Math::floor((p.pos(1)-plo[1])*dxi[1])),
66  p.idata(ERFParticlesIntIdxAoS::k) ) );
67  iv[0] += domain.smallEnd()[0];
68  iv[1] += domain.smallEnd()[1];
69  return iv;
70  }
71 };
72 
73 template <typename P>
74 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
75 void update_location_idata ( P& a_p,
76  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& a_plo,
77  amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& a_dxi,
78  const amrex::Array4<amrex::Real const>& a_height_arr )
79 {
80  amrex::IntVect iv( int(amrex::Math::floor((a_p.pos(0)-a_plo[0])*a_dxi[0])),
81  int(amrex::Math::floor((a_p.pos(1)-a_plo[1])*a_dxi[1])),
82  a_p.idata(ERFParticlesIntIdxAoS::k) );
83 
84  if (a_height_arr) {
85  amrex::Real lx = (a_p.pos(0)-a_plo[0])*a_dxi[0] - static_cast<amrex::Real>(iv[0]);
86  amrex::Real ly = (a_p.pos(1)-a_plo[1])*a_dxi[1] - static_cast<amrex::Real>(iv[1]);
87  auto zlo = a_height_arr(iv[0] ,iv[1] ,iv[2] ) * (1.0-lx) * (1.0-ly) +
88  a_height_arr(iv[0]+1,iv[1] ,iv[2] ) * lx * (1.0-ly) +
89  a_height_arr(iv[0] ,iv[1]+1,iv[2] ) * (1.0-lx) * ly +
90  a_height_arr(iv[0]+1,iv[1]+1,iv[2] ) * lx * ly;
91  auto zhi = a_height_arr(iv[0] ,iv[1] ,iv[2]+1) * (1.0-lx) * (1.0-ly) +
92  a_height_arr(iv[0]+1,iv[1] ,iv[2]+1) * lx * (1.0-ly) +
93  a_height_arr(iv[0] ,iv[1]+1,iv[2]+1) * (1.0-lx) * ly +
94  a_height_arr(iv[0]+1,iv[1]+1,iv[2]+1) * lx * ly;
95 
96  if (a_p.pos(2) > zhi) {
97  a_p.idata(ERFParticlesIntIdxAoS::k) += 1;
98  } else if (a_p.pos(2) <= zlo) {
99  a_p.idata(ERFParticlesIntIdxAoS::k) -= 1;
100  }
101  }
102 }
103 
104 class ERFPC : public amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, // AoS real attributes
105  ERFParticlesIntIdxAoS::ncomps, // AoS integer attributes
106  ERFParticlesRealIdxSoA::ncomps, // SoA real attributes
107  ERFParticlesIntIdxSoA::ncomps, // SoA integer attributes
108  amrex::DefaultAllocator,
109  ERFParticlesAssignor >
110 {
111  public:
112 
113  /*! Constructor */
114  ERFPC ( amrex::ParGDBBase* a_gdb,
115  const std::string& a_name = "particles" )
116  : amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, // AoS real attributes
117  ERFParticlesIntIdxAoS::ncomps, // AoS integer attributes
118  ERFParticlesRealIdxSoA::ncomps, // SoA real attributes
119  ERFParticlesIntIdxSoA::ncomps, // SoA integer attributes
120  amrex::DefaultAllocator,
121  ERFParticlesAssignor> (a_gdb)
122  {
123  BL_PROFILE("ERFPCPC::ERFPC()");
124  m_name = a_name;
125  readInputs();
126  }
127 
128  /*! Constructor */
129  ERFPC ( const amrex::Geometry& a_geom,
130  const amrex::DistributionMapping& a_dmap,
131  const amrex::BoxArray& a_ba,
132  const std::string& a_name = "particles" )
133  : amrex::ParticleContainer< ERFParticlesRealIdxAoS::ncomps, // AoS real attributes
134  ERFParticlesIntIdxAoS::ncomps, // AoS real attributes
135  ERFParticlesRealIdxSoA::ncomps, // SoA real attributes
136  ERFParticlesIntIdxSoA::ncomps, // SoA integer attributes
137  amrex::DefaultAllocator,
138  ERFParticlesAssignor> ( a_geom, a_dmap, a_ba )
139  {
140  BL_PROFILE("ERFPCPC::ERFPC()");
141  m_name = a_name;
142  readInputs();
143  }
144 
145  /*! Initialize particles in domain */
146  virtual void InitializeParticles (const std::unique_ptr<amrex::MultiFab>& a_ptr = nullptr);
147 
148  /*! Evolve particles for one time step */
149  virtual void EvolveParticles ( int,
150  amrex::Real,
151  amrex::Vector<amrex::Vector<amrex::MultiFab>>&,
152  const amrex::Vector<std::unique_ptr<amrex::MultiFab>>& );
153 
154  /*! Get real-type particle attribute names */
155  virtual amrex::Vector<std::string> varNames () const
156  {
157  BL_PROFILE("ERFPCPC::varNames()");
158  return {AMREX_D_DECL("xvel","yvel","zvel"),"mass"};
159  }
160 
161  /*! Get real-type particle attribute names */
162  virtual amrex::Vector<std::string> meshPlotVarNames () const
163  {
164  BL_PROFILE("ERFPCPC::varNames()");
165  return {"mass_density"};
166  }
167 
168  /*! Uses midpoint method to advance particles using flow velocity. */
169  virtual void AdvectWithFlow ( amrex::MultiFab*,
170  int,
171  amrex::Real,
172  const std::unique_ptr<amrex::MultiFab>& );
173 
174  /*! Uses midpoint method to advance particles falling under gravity. */
175  virtual void AdvectWithGravity ( int,
176  amrex::Real,
177  const std::unique_ptr<amrex::MultiFab>& );
178 
179  /*! Compute mass density */
180  virtual void massDensity ( amrex::MultiFab&, const int&, const int& a_comp = 0) const;
181 
182  /*! Compute mesh variable from particles */
183  virtual void computeMeshVar( const std::string& a_var_name,
184  amrex::MultiFab& a_mf,
185  const int a_lev) const
186  {
187  if (a_var_name == "mass_density") {
188  massDensity( a_mf, a_lev );
189  } else {
190  a_mf.setVal(0.0);
191  }
192  }
193 
194  /*! Specify if particles should advect with flow */
195  inline void setAdvectWithFlow (bool a_flag)
196  {
197  BL_PROFILE("ERFPCPC::setAdvectWithFlow()");
198  m_advect_w_flow = a_flag;
199  }
200  /*! Specify if particles fall under gravity */
201  inline void setAdvectWithGravity (bool a_flag)
202  {
203  BL_PROFILE("ERFPCPC::setAdvectWithGravity()");
204  m_advect_w_gravity = a_flag;
205  }
206 
207  // the following functions should ideally be private or protected, but need to be
208  // public due to CUDA extended lambda capture rules
209 
210  /*! Default particle initialization */
211  void initializeParticlesUniformDistributionInBox (const std::unique_ptr<amrex::MultiFab>& a_ptr,
212  const amrex::RealBox& particle_box);
213 
214  protected:
215 
216  bool m_advect_w_flow; /*!< advect with flow velocity */
217  bool m_advect_w_gravity; /*!< advect under gravitational force */
218 
219  amrex::RealBox m_particle_box; /*!< box within which to place particles */
220 
221  std::string m_name; /*!< name of this particle species */
222 
223  std::string m_initialization_type; /*!< initial particle distribution type */
224  int m_ppc_init; /*!< initial number of particles per cell */
225 
226  bool m_stable_redistribute; /*!< use stable redistribute for deterministic simulations */
227 
228  /*! read inputs from file */
229  virtual void readInputs ();
230 
231  private:
232 
233  bool place_randomly_in_cells; /*!< place particles at random positions? */
234 };
235 
236 #endif
237 #endif
amrex
Definition: ERF_ConsoleIO.cpp:12