REMORA
Regional Modeling of Oceans Refined Adaptively
Loading...
Searching...
No Matches
REMORA_PC.H
Go to the documentation of this file.
1#ifndef REMORA_PC_H_
2#define REMORA_PC_H_
3
4#ifdef REMORA_USE_PARTICLES
5
6#include <string>
7#include <AMReX_Particles.H>
8
9struct REMORAParticlesIntIdxAoS
10{
11 enum {
12 k = 0,
13 ncomps
14 };
15};
16
17struct REMORAParticlesRealIdxAoS
18{
19 enum {
20 ncomps = 0
21 };
22};
23
24struct REMORAParticlesIntIdxSoA
25{
26 enum {
27 ncomps = 0
28 };
29};
30
31struct REMORAParticlesRealIdxSoA
32{
33 enum {
34 vx = 0,
35 vy,
36 vz,
37 mass,
38 ncomps
39 };
40};
41
42namespace REMORAParticleInitializations
43{
44 /* list of particle initializations */
45 const std::string init_box_uniform = "box";
46}
47
48namespace REMORAParticleNames
49{
50 const std::string tracers = "tracer_particles";
51 const std::string hydro = "hydro_particles";
52}
53
54struct REMORAParticlesAssignor
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(REMORAParticlesIntIdxAoS::k) ) );
67 iv[0] += domain.smallEnd()[0];
68 iv[1] += domain.smallEnd()[1];
69 return iv;
70 }
71};
72
73template <typename P>
74AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
75void 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(REMORAParticlesIntIdxAoS::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(REMORAParticlesIntIdxAoS::k) += 1;
98 } else if (a_p.pos(2) <= zlo) {
99 a_p.idata(REMORAParticlesIntIdxAoS::k) -= 1;
100 }
101 }
102}
103
104class REMORAPC : public amrex::ParticleContainer< REMORAParticlesRealIdxAoS::ncomps, // AoS real attributes
105 REMORAParticlesIntIdxAoS::ncomps, // AoS integer attributes
106 REMORAParticlesRealIdxSoA::ncomps, // SoA real attributes
107 REMORAParticlesIntIdxSoA::ncomps, // SoA integer attributes
108 amrex::DefaultAllocator,
109 REMORAParticlesAssignor >
110{
111 public:
112
113 /*! Constructor */
114 REMORAPC ( amrex::ParGDBBase* a_gdb,
115 const std::string& a_name = "particles" )
116 : amrex::ParticleContainer< REMORAParticlesRealIdxAoS::ncomps, // AoS real attributes
117 REMORAParticlesIntIdxAoS::ncomps, // AoS integer attributes
118 REMORAParticlesRealIdxSoA::ncomps, // SoA real attributes
119 REMORAParticlesIntIdxSoA::ncomps, // SoA integer attributes
120 amrex::DefaultAllocator,
121 REMORAParticlesAssignor> (a_gdb)
122 {
123 BL_PROFILE("REMORAPCPC::REMORAPC()");
124 m_name = a_name;
125 readInputs();
126 }
127
128 /*! Constructor */
129 REMORAPC ( 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< REMORAParticlesRealIdxAoS::ncomps, // AoS real attributes
134 REMORAParticlesIntIdxAoS::ncomps, // AoS real attributes
135 REMORAParticlesRealIdxSoA::ncomps, // SoA real attributes
136 REMORAParticlesIntIdxSoA::ncomps, // SoA integer attributes
137 amrex::DefaultAllocator,
138 REMORAParticlesAssignor> ( a_geom, a_dmap, a_ba )
139 {
140 BL_PROFILE("REMORAPCPC::REMORAPC()");
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 lev,
150 amrex::Real dt,
151 amrex::Vector<amrex::MultiFab const*>& a_flow_vel,
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("REMORAPCPC::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("REMORAPCPC::varNames()");
165 return {"mass_density"};
166 }
167
168 /*! Uses midpoint method to advance particles using flow velocity. */
169 virtual void AdvectWithFlow ( int lev,
170 amrex::Real dt,
171 amrex::Vector<amrex::MultiFab const*>&,
172 const std::unique_ptr<amrex::MultiFab>& );
173
174 /*! Compute mass density */
175 virtual void massDensity ( amrex::MultiFab&, const int&, const int& a_comp = 0) const;
176
177 /*! Compute mesh variable from particles */
178 virtual void computeMeshVar( const std::string& a_var_name,
179 amrex::MultiFab& a_mf,
180 const int a_lev) const
181 {
182 if (a_var_name == "mass_density") {
183 massDensity( a_mf, a_lev );
184 } else {
185 a_mf.setVal(0.0);
186 }
187 }
188
189 /*! Specify if particles should advect with flow */
190 inline void setAdvectWithFlow (bool a_flag)
191 {
192 BL_PROFILE("REMORAPCPC::setAdvectWithFlow()");
193 m_advect_w_flow = a_flag;
194 }
195
196 // the following functions should ideally be private or protected, but need to be
197 // public due to CUDA extended lambda capture rules
198
199 /*! Default particle initialization */
200 void initializeParticlesUniformDistributionInBox (const std::unique_ptr<amrex::MultiFab>& a_ptr,
201 const amrex::RealBox& particle_box);
202
203 protected:
204
205 bool m_advect_w_flow; /*!< advect with flow velocity */
206
207 amrex::RealBox m_particle_box; /*!< box within which to place particles */
208
209 std::string m_name; /*!< name of this particle species */
210
211 std::string m_initialization_type; /*!< initial particle distribution type */
212 int m_ppc_init; /*!< initial number of particles per cell */
213
214 /*! read inputs from file */
215 virtual void readInputs ();
216
217 private:
218
219 bool place_randomly_in_cells; /*!< place particles at random positions? */
220};
221
222#endif
223#endif