Block-Structured AMR Software Framework
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
AMReX_AmrParticles.H
Go to the documentation of this file.
1#ifndef AMREX_AmrParticles_H_
2#define AMREX_AmrParticles_H_
3#include <AMReX_Config.H>
4
5#include <AMReX_Particles.H>
7#include <AMReX_AmrParGDB.H>
10
11namespace amrex {
12
13template <typename ParticleType, int NArrayReal, int NArrayInt,
14 template<class> class Allocator, class CellAssignor>
15void
16ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
17::AssignDensity (int rho_index,
18 Vector<std::unique_ptr<MultiFab> >& mf_to_be_filled,
19 int lev_min, int ncomp, int finest_level, int ngrow) const
20{
21
22 BL_PROFILE("ParticleContainer::AssignDensity()");
23
24 if (rho_index != 0) { amrex::Abort("AssignDensity only works if rho_index = 0"); }
25
26 BL_ASSERT(NStructReal >= 1);
27 BL_ASSERT(NStructReal >= ncomp);
28 BL_ASSERT(ncomp == 1 || ncomp == AMREX_SPACEDIM+1);
29
30 if (finest_level == -1) {
31 finest_level = this->finestLevel();
32 }
33 while (!this->m_gdb->LevelDefined(finest_level)) {
34 finest_level--;
35 }
36
37 ngrow = std::max(ngrow, 2);
38
39 // Create the space for mf_to_be_filled, regardless of whether we'll need a temporary mf
40 mf_to_be_filled.resize(finest_level+1);
41 for (int lev = lev_min; lev <= finest_level; lev++)
42 {
43 auto ng = lev == lev_min ? IntVect(AMREX_D_DECL(ngrow,ngrow,ngrow)) : this->m_gdb->refRatio(lev-1);
44 mf_to_be_filled[lev] = std::make_unique<MultiFab>(this->m_gdb->boxArray(lev),
45 this->m_gdb->DistributionMap(lev),
46 ncomp, ng);
47 mf_to_be_filled[lev]->setVal(0.0);
48 }
49
50 // Test whether the grid structure of the boxArray is the same as the ParticleBoxArray at all levels
51 bool all_grids_the_same = true;
52 for (int lev = lev_min; lev <= finest_level; lev++)
53 {
54 if (!this->OnSameGrids(lev, *mf_to_be_filled[lev]))
55 {
56 all_grids_the_same = false;
57 break;
58 }
59 }
60
62 if (!all_grids_the_same)
63 {
64 // Create the space for the temporary, mf_part
65 mf_part.resize(finest_level+1);
66 for (int lev = lev_min; lev <= finest_level; lev++)
67 {
68 auto ng = lev == lev_min ? IntVect(AMREX_D_DECL(ngrow,ngrow,ngrow)) : this->m_gdb->refRatio(lev-1);
69 mf_part[lev] = std::make_unique<MultiFab>(this->ParticleBoxArray(lev),
70 this->ParticleDistributionMap(lev),
71 ncomp, ng);
72 mf_part[lev]->setVal(0.0);
73 }
74 }
75
76 auto & mf = (all_grids_the_same) ? mf_to_be_filled : mf_part;
77
78 if (finest_level == 0)
79 {
80 //
81 // Just use the far simpler single-level version.
82 //
83 this->AssignCellDensitySingleLevel(rho_index, *mf[0], 0, ncomp);
84 if ( ! all_grids_the_same) {
85 mf_to_be_filled[0]->ParallelCopy(*mf[0],0,0,ncomp,0,0);
86 }
87 return;
88 }
89
90 // configure this to do a no-op at the physical boundaries.
91 int lo_bc[] = {BCType::int_dir, BCType::int_dir, BCType::int_dir};
92 int hi_bc[] = {BCType::int_dir, BCType::int_dir, BCType::int_dir};
93 Vector<BCRec> bcs(ncomp, BCRec(lo_bc, hi_bc));
94 PCInterp mapper;
95
96 Vector<std::unique_ptr<MultiFab> > tmp(finest_level+1);
97 for (int lev = lev_min; lev <= finest_level; ++lev) {
98 const BoxArray& ba = mf[lev]->boxArray();
99 const DistributionMapping& dm = mf[lev]->DistributionMap();
100 tmp[lev] = std::make_unique<MultiFab>(ba, dm, ncomp, 0);
101 tmp[lev]->setVal(0.0);
102 }
103
104 for (int lev = lev_min; lev <= finest_level; ++lev) {
105
106 this->AssignCellDensitySingleLevel(rho_index, *mf[lev], lev, ncomp, 0);
107
108 if (lev < finest_level) {
109 PhysBCFunctNoOp cphysbc, fphysbc;
110 amrex::InterpFromCoarseLevel(*tmp[lev+1], 0.0, *mf[lev],
111 rho_index, rho_index, ncomp,
112 this->m_gdb->Geom(lev),
113 this->m_gdb->Geom(lev+1),
114 cphysbc, 0, fphysbc, 0,
115 this->m_gdb->refRatio(lev), &mapper,
116 bcs, 0);
117 }
118
119 if (lev > lev_min) {
120 // Note - this will double count the mass on the coarse level in
121 // regions covered by the fine level, but this will be corrected
122 // below in the call to average_down.
123 amrex::sum_fine_to_coarse(*mf[lev], *mf[lev-1], rho_index, ncomp,
124 this->m_gdb->refRatio(lev-1),
125 this->m_gdb->Geom(lev-1),
126 this->m_gdb->Geom(lev));
127 }
128
129 mf[lev]->plus(*tmp[lev], rho_index, ncomp, 0);
130 }
131
132 for (int lev = finest_level - 1; lev >= lev_min; --lev) {
133 amrex::average_down(*mf[lev+1], *mf[lev], rho_index, ncomp, this->m_gdb->refRatio(lev));
134 }
135
136 if (!all_grids_the_same) {
137 for (int lev = lev_min; lev <= finest_level; lev++) {
138 //
139 // We haven't interpolated the ghost cells so we can't copy them
140 //
141 mf_to_be_filled[lev]->ParallelCopy(*mf_part[lev],0,0,ncomp,0,0);
142 }
143 }
144 if (lev_min > 0) {
145 int nlevels = finest_level - lev_min + 1;
146 for (int i = 0; i < nlevels; i++)
147 {
148 mf_to_be_filled[i] = std::move(mf_to_be_filled[i+lev_min]);
149 }
150 mf_to_be_filled.resize(nlevels);
151 }
152}
153
154template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
155void
156ParticleToMesh (PC const& pc, const Vector<MultiFab*>& mf,
157 int lev_min, int lev_max, F&& f,
158 bool zero_out_input=true, bool vol_weight=true)
159{
160 BL_PROFILE("amrex::ParticleToMesh");
161
162 if (lev_max == -1) { lev_max = pc.finestLevel(); }
163 while (!pc.GetParGDB()->LevelDefined(lev_max)) { lev_max--; }
164
165 if (lev_max == 0)
166 {
167 ParticleToMesh(pc, *mf[0], 0, std::forward<F>(f), zero_out_input);
168 if (vol_weight) {
169 const Real* dx = pc.Geom(0).CellSize();
170 const Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]);
171 mf[0]->mult(Real(1.0)/vol, 0, mf[0]->nComp(), mf[0]->nGrow());
172 }
173 return;
174 }
175
176 int ngrow = amrex::max(AMREX_D_DECL(mf[0]->nGrow(0),
177 mf[0]->nGrow(1),
178 mf[0]->nGrow(2)), 2);
179
180 if (zero_out_input)
181 {
182 for (int lev = lev_min; lev <= lev_max; ++lev)
183 {
184 mf[lev]->setVal(Real(0.0));
185 }
186 }
187
188 Vector<MultiFab> mf_part(lev_max+1);
189 Vector<MultiFab> mf_tmp( lev_max+1);
190 for (int lev = lev_min; lev <= lev_max; ++lev)
191 {
192 mf_part[lev].define(pc.ParticleBoxArray(lev),
193 pc.ParticleDistributionMap(lev),
194 mf[lev]->nComp(), ngrow);
195 mf_tmp[lev].define( pc.ParticleBoxArray(lev),
196 pc.ParticleDistributionMap(lev),
197 mf[lev]->nComp(), 0);
198 mf_part[lev].setVal(0.0);
199 mf_tmp[lev].setVal( 0.0);
200 }
201
202 // configure this to do a no-op at the physical boundaries.
203 int lo_bc[] = {BCType::int_dir, BCType::int_dir, BCType::int_dir};
204 int hi_bc[] = {BCType::int_dir, BCType::int_dir, BCType::int_dir};
205 Vector<BCRec> bcs(mf_part[lev_min].nComp(), BCRec(lo_bc, hi_bc));
206 PCInterp mapper;
207
208 for (int lev = lev_min; lev <= lev_max; ++lev)
209 {
210 ParticleToMesh(pc, mf_part[lev], lev, f, zero_out_input);
211 if (vol_weight) {
212 const Real* dx = pc.Geom(lev).CellSize();
213 const Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]);
214 mf_part[lev].mult(Real(1.0)/vol, 0, mf_part[lev].nComp(), mf_part[lev].nGrow());
215 }
216
217 if (lev < lev_max) {
218 PhysBCFunctNoOp cphysbc, fphysbc;
219 amrex::InterpFromCoarseLevel(mf_tmp[lev+1], 0.0, mf_part[lev],
220 0, 0, mf_part[lev].nComp(),
221 pc.GetParGDB()->Geom(lev),
222 pc.GetParGDB()->Geom(lev+1),
223 cphysbc, 0, fphysbc, 0,
224 pc.GetParGDB()->refRatio(lev), &mapper,
225 bcs, 0);
226 }
227
228 if (lev > lev_min) {
229 // Note - this will double count the mass on the coarse level in
230 // regions covered by the fine level, but this will be corrected
231 // below in the call to average_down.
232 amrex::sum_fine_to_coarse(mf_part[lev], mf_part[lev-1], 0, mf_part[lev].nComp(),
233 pc.GetParGDB()->refRatio(lev-1),
234 pc.GetParGDB()->Geom(lev-1),
235 pc.GetParGDB()->Geom(lev));
236 }
237
238 mf_part[lev].plus(mf_tmp[lev], 0, mf_part[lev].nComp(), 0);
239 }
240
241 for (int lev = lev_max-1; lev >= lev_min; --lev) {
242 amrex::average_down(mf_part[lev+1], mf_part[lev], 0, mf_part[lev].nComp(),
243 pc.GetParGDB()->refRatio(lev));
244 }
245
246 for (int lev = lev_min; lev <= lev_max; lev++)
247 {
248 mf[lev]->ParallelCopy(mf_part[lev],0,0,mf_part[lev].nComp(),0,0);
249 }
250}
251
252template <typename T_ParticleType, int NArrayReal=0, int NArrayInt=0,
253 template<class> class Allocator=DefaultAllocator, class CellAssignor=DefaultAssignor>
254class AmrParticleContainer_impl // NOLINT(cppcoreguidelines-virtual-class-destructor)
255 : public ParticleContainer_impl<T_ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
256{
257
258public:
259
260 using ParticleType = T_ParticleType;
261
266
271
273 const Vector<DistributionMapping> & dmap,
274 const Vector<BoxArray> & ba,
275 const Vector<int> & rr)
276 : ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>(geom, dmap, ba, rr)
277 {
278 }
279
280 ~AmrParticleContainer_impl () override = default;
281
284
286 AmrParticleContainer_impl& operator= ( AmrParticleContainer_impl && ) noexcept = default;
287};
288
289template <int T_NStructReal, int T_NStructInt=0, int T_NArrayReal=0, int T_NArrayInt=0,
290 template<class> class Allocator=DefaultAllocator, class CellAssignor=DefaultAssignor>
291using AmrParticleContainer = AmrParticleContainer_impl<Particle<T_NStructReal, T_NStructInt>, T_NArrayReal, T_NArrayInt, Allocator, CellAssignor>;
292
311
312}
313
314#endif
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define BL_ASSERT(EX)
Definition AMReX_BLassert.H:39
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:129
Provide basic functionalities to set up an AMR hierarchy.
Definition AMReX_AmrCore.H:25
Definition AMReX_AmrParticles.H:256
AmrParticleContainer_impl(const AmrParticleContainer_impl &)=delete
AmrParticleContainer_impl & operator=(const AmrParticleContainer_impl &)=delete
AmrParticleContainer_impl(const Vector< Geometry > &geom, const Vector< DistributionMapping > &dmap, const Vector< BoxArray > &ba, const Vector< int > &rr)
Definition AMReX_AmrParticles.H:272
~AmrParticleContainer_impl() override=default
AmrParticleContainer_impl(AmrCore *amr_core)
Definition AMReX_AmrParticles.H:267
AmrParticleContainer_impl(AmrParticleContainer_impl &&) noexcept=default
T_ParticleType ParticleType
Definition AMReX_AmrParticles.H:260
AmrParticleContainer_impl()
Definition AMReX_AmrParticles.H:262
Definition AMReX_AmrParticles.H:295
AmrTracerParticleContainer(AmrTracerParticleContainer &&) noexcept=default
~AmrTracerParticleContainer() override=default
AmrTracerParticleContainer(const AmrTracerParticleContainer &)=delete
AmrTracerParticleContainer(AmrCore *amr_core)
Definition AMReX_AmrParticles.H:298
Definition AMReX_GpuAllocators.H:114
Boundary Condition Records. Necessary information and functions for computing boundary conditions.
Definition AMReX_BCRec.H:17
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:550
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:41
Piecewise Constant interpolation on cell centered data.
Definition AMReX_Interpolater.H:451
const ParGDBBase * GetParGDB() const
Get the ParGDB object used to define this container (const version)
Definition AMReX_ParticleContainerBase.H:225
A distributed container for Particles sorted onto the levels, grids, and tiles of a block-structured ...
Definition AMReX_ParticleContainer.H:146
static constexpr int NArrayInt
Number of extra integer components stored in struct-of-array form.
Definition AMReX_ParticleContainer.H:159
static constexpr int NArrayReal
Number of extra Real components stored in struct-of-array form.
Definition AMReX_ParticleContainer.H:157
T_CellAssignor CellAssignor
Definition AMReX_ParticleContainer.H:150
Definition AMReX_PhysBCFunct.H:120
Definition AMReX_TracerParticles.H:11
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:27
Definition AMReX_Amr.cpp:49
int nComp(FabArrayBase const &fa)
amrex::ArenaAllocator< T > DefaultAllocator
Definition AMReX_GpuAllocators.H:194
void average_down(const MultiFab &S_fine, MultiFab &S_crse, const Geometry &fgeom, const Geometry &cgeom, int scomp, int ncomp, int rr)
Definition AMReX_MultiFabUtil.cpp:314
std::enable_if_t< IsFabArray< MF >::value > InterpFromCoarseLevel(MF &mf, Real time, const MF &cmf, int scomp, int dcomp, int ncomp, const Geometry &cgeom, const Geometry &fgeom, BC &cbc, int cbccomp, BC &fbc, int fbccomp, const IntVect &ratio, Interp *mapper, const Vector< BCRec > &bcs, int bcscomp, const PreInterpHook &pre_interp={}, const PostInterpHook &post_interp={})
Fill with interpolation of coarse level data.
Definition AMReX_FillPatchUtil_I.H:984
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
void ParticleToMesh(PC const &pc, const Vector< MultiFab * > &mf, int lev_min, int lev_max, F &&f, bool zero_out_input=true, bool vol_weight=true)
Definition AMReX_AmrParticles.H:156
void sum_fine_to_coarse(const MultiFab &S_fine, MultiFab &S_crse, int scomp, int ncomp, const IntVect &ratio, const Geometry &cgeom, const Geometry &)
Definition AMReX_MultiFabUtil.cpp:393
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE constexpr const T & max(const T &a, const T &b) noexcept
Definition AMReX_Algorithm.H:35
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:230
Definition AMReX_ParticleUtil.H:432
The struct used to store particles.
Definition AMReX_Particle.H:295