Block-Structured AMR Software Framework
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>
8 #include <AMReX_Interpolater.H>
9 #include <AMReX_FillPatchUtil.H>
10 
11 namespace amrex {
12 
13 template <typename ParticleType, int NArrayReal, int NArrayInt,
14  template<class> class Allocator, class CellAssignor>
15 void
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 
154 template <class PC, class F, std::enable_if_t<IsParticleContainer<PC>::value, int> foo = 0>
155 void
156 ParticleToMesh (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 
252 template <typename T_ParticleType, int NArrayReal=0, int NArrayInt=0,
253  template<class> class Allocator=DefaultAllocator, class CellAssignor=DefaultAssignor>
254 class AmrParticleContainer_impl // NOLINT(cppcoreguidelines-virtual-class-destructor)
255  : public ParticleContainer_impl<T_ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
256 {
257 
258 public:
259 
260  using ParticleType = T_ParticleType;
261 
264  {
265  }
266 
269  {
270  }
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 
289 template <int T_NStructReal, int T_NStructInt=0, int T_NArrayReal=0, int T_NArrayInt=0,
290  template<class> class Allocator=DefaultAllocator, class CellAssignor=DefaultAssignor>
291 using AmrParticleContainer = AmrParticleContainer_impl<Particle<T_NStructReal, T_NStructInt>, T_NArrayReal, T_NArrayInt, Allocator, CellAssignor>;
292 
294  : public TracerParticleContainer
295 {
296 public:
297 
299  : TracerParticleContainer(amr_core->GetParGDB())
300  {
301  }
302 
303  ~AmrTracerParticleContainer () override = default;
304 
307 
309  AmrTracerParticleContainer& operator= ( AmrTracerParticleContainer && ) noexcept = default;
310 };
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
#define AMREX_D_DECL(a, b, c)
Definition: AMReX_SPACE.H:104
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:145
void AssignDensity(int rho_index, Vector< std::unique_ptr< MultiFab > > &mf_to_be_filled, int lev_min, int ncomp, int finest_level, int ngrow=2) const
Functions depending the layout of the data. Use with caution.
Definition: AMReX_AmrParticles.H:17
static constexpr int NArrayInt
Number of extra integer components stored in struct-of-array form.
Definition: AMReX_ParticleContainer.H:158
static constexpr int NArrayReal
Number of extra Real components stored in struct-of-array form.
Definition: AMReX_ParticleContainer.H:156
DefaultAssignor CellAssignor
Definition: AMReX_ParticleContainer.H:149
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
static int f(amrex::Real t, N_Vector y_data, N_Vector y_rhs, void *user_data)
Definition: AMReX_SundialsIntegrator.H:44
@ max
Definition: AMReX_ParallelReduce.H:17
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
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
Definition: AMReX_Algorithm.H:35
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
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition: AMReX.cpp:225
Definition: AMReX_ParticleUtil.H:432
The struct used to store particles.
Definition: AMReX_Particle.H:295