1 #ifndef BOXLIB_PARTICLE_HPP
2 #define BOXLIB_PARTICLE_HPP
6 #include <AMReX_BLFort.H>
7 #include <AMReX_MultiFabUtil.H>
8 #include <AMReX_MultiFabUtil_F.H>
9 #include <AMReX_Interpolater.H>
10 #include <AMReX_FillPatchUtil.H>
12 template<
class PLayout>
19 template<
class PLayout>
26 template<
class PLayout>
27 template<
class FT,
unsigned Dim,
class PT>
33 if ( lbase == lfine ) {
34 this->
scatter(attrib, *(f[lbase].
get()), pp, pbin, bin, lbase);
38 const PLayout *layout_p = &this->getLayout();
39 int nGrow = layout_p->refRatio(lbase)[0];
42 for (
int lev = lbase; lev <= lfine; ++lev) {
44 f[lev]->setVal(0.0, f[lev]->nGrow());
48 tmp[lev].reset(
new AmrField_t(ba, dm, 1, nGrow));
49 tmp[lev]->setVal(0.0, nGrow);
52 this->AssignDensityFort(attrib, tmp, lbase, 1, lfine, pbin, bin);
54 for (
int lev = lbase; lev <= lfine; ++lev)
55 AmrField_t::Copy(*f[lev], *tmp[lev], 0, 0, f[lev]->nComp(), f[lev]->nGrow());
59 template<
class PLayout>
60 template <
class FT,
unsigned Dim,
class PT>
72 this->AssignCellDensitySingleLevelFort(attrib, tmp, level, pbin, bin);
74 f.setVal(0.0, f.nGrow());
76 AmrField_t::Copy(f, tmp, 0, 0, f.nComp(), f.nGrow());
80 template<
class PLayout>
81 template<
class FT,
unsigned Dim,
class PT>
86 this->InterpolateFort(attrib, f, lbase, lfine);
92 template<
class PLayout>
93 template <
class AType>
96 int lev_min,
int ncomp,
int finest_level,
101 const PLayout *layout_p = &this->getLayout();
106 amrex::PhysBCFunct cphysbc, fphysbc;
107 int lo_bc[] = {INT_DIR, INT_DIR, INT_DIR};
108 int hi_bc[] = {INT_DIR, INT_DIR, INT_DIR};
109 amrex::Vector<amrex::BCRec> bcs(1, amrex::BCRec(lo_bc, hi_bc));
110 amrex::PCInterp mapper;
113 for (
int lev = lev_min; lev <= finest_level; ++lev) {
114 const AmrGrid_t& ba = mf_to_be_filled[lev]->boxArray();
115 const AmrProcMap_t& dm = mf_to_be_filled[lev]->DistributionMap();
117 tmp[lev]->setVal(0.0);
120 for (
int lev = lev_min; lev <= finest_level; ++lev) {
121 AssignCellDensitySingleLevelFort(pa, *mf_to_be_filled[lev], lev, pbin, bin, 1, 0);
123 if (lev < finest_level) {
124 amrex::InterpFromCoarseLevel(*tmp[lev+1], 0.0, *mf_to_be_filled[lev],
125 rho_index, rho_index, ncomp,
126 layout_p->Geom(lev), layout_p->Geom(lev+1),
128 layout_p->refRatio(lev), &mapper, bcs);
135 amrex::sum_fine_to_coarse(*mf_to_be_filled[lev],
136 *mf_to_be_filled[lev-1],
137 rho_index, 1, layout_p->refRatio(lev-1),
138 layout_p->Geom(lev-1), layout_p->Geom(lev));
141 mf_to_be_filled[lev]->plus(*tmp[lev], rho_index, ncomp, 0);
144 for (
int lev = finest_level - 1; lev >= lev_min; --lev) {
145 amrex::average_down(*mf_to_be_filled[lev+1],
146 *mf_to_be_filled[lev], rho_index, ncomp, layout_p->refRatio(lev));
154 template<
class PLayout>
155 template <
class AType>
162 int particle_lvl_offset)
const
166 const PLayout *layout_p = &this->getLayout();
170 if (layout_p->OnSameGrids(level, mf_to_be_filled)) {
173 mf_pointer = &mf_to_be_filled;
178 mf_pointer =
new AmrField_t(layout_p->ParticleBoxArray(level),
179 layout_p->ParticleDistributionMap(level),
180 ncomp, mf_to_be_filled.nGrow());
187 if (mf_pointer->nGrow() < 1)
188 throw OpalException(
"BoxLibParticle::AssignCellDensitySingleLevelFort()",
189 "Must have at least one ghost cell when in AssignCellDensitySingleLevelFort");
196 if (gm.isAnyPeriodic() && ! gm.isAllPeriodic()) {
197 throw OpalException(
"BoxLibParticle::AssignCellDensitySingleLevelFort()",
198 "Problem must be periodic in no or all directions");
201 for (amrex::MFIter mfi(*mf_pointer); mfi.isValid(); ++mfi) {
202 (*mf_pointer)[mfi].setVal(0);
207 size_t lBegin = LocalNumPerLevel.
begin(level);
208 size_t lEnd = LocalNumPerLevel.
end(level);
210 AmrReal_t inv_dx[3] = { 1.0 / dx[0], 1.0 / dx[1], 1.0 / dx[2] };
211 double lxyz[3] = { 0.0, 0.0, 0.0 };
212 double wxyz_hi[3] = { 0.0, 0.0, 0.0 };
213 double wxyz_lo[3] = { 0.0, 0.0, 0.0 };
214 int ijk[3] = {0, 0, 0};
216 for (
size_t ip = lBegin; ip < lEnd; ++ip) {
218 if ( bin > -1 && pbin[ip] != bin )
221 const int grid = this->Grid[ip];
226 for (
int i = 0; i < 3; ++i) {
227 lxyz[i] = ( this->
R[ip](i) - plo[i] ) * inv_dx[i] + 0.5;
229 wxyz_hi[i] = lxyz[i] - ijk[i];
230 wxyz_lo[i] = 1.0 - wxyz_hi[i];
246 fab(i1, 0) += wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*pa[ip];
247 fab(i2, 0) += wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*pa[ip];
248 fab(i3, 0) += wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*pa[ip];
249 fab(i4, 0) += wxyz_lo[0]*wxyz_hi[1]*wxyz_hi[2]*pa[ip];
250 fab(i5, 0) += wxyz_hi[0]*wxyz_lo[1]*wxyz_lo[2]*pa[ip];
251 fab(i6, 0) += wxyz_hi[0]*wxyz_lo[1]*wxyz_hi[2]*pa[ip];
252 fab(i7, 0) += wxyz_hi[0]*wxyz_hi[1]*wxyz_lo[2]*pa[ip];
253 fab(i8, 0) += wxyz_hi[0]*wxyz_hi[1]*wxyz_hi[2]*pa[ip];
257 mf_pointer->SumBoundary(gm.periodicity());
262 const AmrReal_t vol = D_TERM(dx[0], *dx[1], *dx[2]);
264 mf_pointer->mult(1.0/vol, 0, 1, mf_pointer->nGrow());
269 if (mf_pointer != &mf_to_be_filled) {
270 mf_to_be_filled.copy(*mf_pointer,0,0,ncomp);
276 template<
class PLayout>
277 template <
class AType>
280 int lev_min,
int lev_max)
282 for (
int lev = lev_min; lev <= lev_max; ++lev) {
284 InterpolateMultiLevelFort(pa, mesh_data, lev);
286 InterpolateSingleLevelFort(pa, mesh_data[lev], lev);
291 template<
class PLayout>
292 template <
class AType>
296 for (std::size_t i = 0; i < mesh_data.size(); ++i) {
297 if (mesh_data[i]->nGrow() < 1)
298 throw OpalException(
"BoxLibParticle::InterpolateSingleLevelFort()",
299 "Must have at least one ghost cell when in InterpolateSingleLevelFort");
302 PLayout *layout_p = &this->getLayout();
310 size_t lBegin = LocalNumPerLevel.
begin(lev);
311 size_t lEnd = LocalNumPerLevel.
end(lev);
314 for (std::size_t i = 0; i < mesh_data.size(); ++i) {
315 mesh_data[i]->FillBoundary(gm.periodicity());
318 AmrReal_t inv_dx[3] = { 1.0 / dx[0], 1.0 / dx[1], 1.0 / dx[2] };
319 double lxyz[3] = { 0.0, 0.0, 0.0 };
320 double wxyz_hi[3] = { 0.0, 0.0, 0.0 };
321 double wxyz_lo[3] = { 0.0, 0.0, 0.0 };
322 int ijk[3] = {0, 0, 0};
323 for (
size_t ip = lBegin; ip < lEnd; ++ip) {
325 const int grid = this->Grid[ip];
333 for (
int i = 0; i < 3; ++i) {
334 lxyz[i] = ( this->
R[ip](i) - plo[i] ) * inv_dx[i] + 0.5;
336 wxyz_hi[i] = lxyz[i] - ijk[i];
337 wxyz_lo[i] = 1.0 - wxyz_hi[i];
353 pa[ip](0) = wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*exfab(i1) +
354 wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*exfab(i2) +
355 wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*exfab(i3) +
356 wxyz_lo[0]*wxyz_hi[1]*wxyz_hi[2]*exfab(i4) +
357 wxyz_hi[0]*wxyz_lo[1]*wxyz_lo[2]*exfab(i5) +
358 wxyz_hi[0]*wxyz_lo[1]*wxyz_hi[2]*exfab(i6) +
359 wxyz_hi[0]*wxyz_hi[1]*wxyz_lo[2]*exfab(i7) +
360 wxyz_hi[0]*wxyz_hi[1]*wxyz_hi[2]*exfab(i8);
362 pa[ip](1) = wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*eyfab(i1) +
363 wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*eyfab(i2) +
364 wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*eyfab(i3) +
365 wxyz_lo[0]*wxyz_hi[1]*wxyz_hi[2]*eyfab(i4) +
366 wxyz_hi[0]*wxyz_lo[1]*wxyz_lo[2]*eyfab(i5) +
367 wxyz_hi[0]*wxyz_lo[1]*wxyz_hi[2]*eyfab(i6) +
368 wxyz_hi[0]*wxyz_hi[1]*wxyz_lo[2]*eyfab(i7) +
369 wxyz_hi[0]*wxyz_hi[1]*wxyz_hi[2]*eyfab(i8);
371 pa[ip](2) = wxyz_lo[0]*wxyz_lo[1]*wxyz_lo[2]*ezfab(i1) +
372 wxyz_lo[0]*wxyz_lo[1]*wxyz_hi[2]*ezfab(i2) +
373 wxyz_lo[0]*wxyz_hi[1]*wxyz_lo[2]*ezfab(i3) +
374 wxyz_lo[0]*wxyz_hi[1]*wxyz_hi[2]*ezfab(i4) +
375 wxyz_hi[0]*wxyz_lo[1]*wxyz_lo[2]*ezfab(i5) +
376 wxyz_hi[0]*wxyz_lo[1]*wxyz_hi[2]*ezfab(i6) +
377 wxyz_hi[0]*wxyz_hi[1]*wxyz_lo[2]*ezfab(i7) +
378 wxyz_hi[0]*wxyz_hi[1]*wxyz_hi[2]*ezfab(i8);
384 template<
class PLayout>
385 template <
class AType>
390 for (std::size_t i = 0; i < mesh_data[lev].size(); ++i) {
391 if (mesh_data[lev][i]->nGrow() < 1)
392 throw OpalException(
"BoxLibParticle::InterpolateMultiLevelFort()",
393 "Must have at least one ghost cell when in InterpolateMultiLevelFort");
396 PLayout *layout_p = &this->getLayout();
401 const AmrReal_t* cdx = layout_p->Geom(lev-1).CellSize();
402 const AmrReal_t* cplo = layout_p->Geom(lev-1).ProbLo();
406 size_t lBegin = LocalNumPerLevel.
begin(lev);
407 size_t lEnd = LocalNumPerLevel.
end(lev);
410 for (std::size_t i = 0; i < mesh_data[lev].size(); ++i) {
411 mesh_data[lev][i]->FillBoundary(gm.periodicity());
414 AmrReal_t inv_fdx[3] = { 1.0 / fdx[0], 1.0 / fdx[1], 1.0 / fdx[2] };
415 AmrReal_t inv_cdx[3] = { 1.0 / cdx[0], 1.0 / cdx[1], 1.0 / cdx[2] };
416 double lxyz[3] = { 0.0, 0.0, 0.0 };
417 double wxyz_hi[3] = { 0.0, 0.0, 0.0 };
418 double wxyz_lo[3] = { 0.0, 0.0, 0.0 };
419 int ijk[3] = { 0, 0, 0 };
421 const AmrGrid_t& fba = mesh_data[lev][0]->boxArray();
422 const AmrProcMap_t& fdmap = mesh_data[lev][0]->DistributionMap();
427 mesh_data[lev][0]->nComp(),
428 mesh_data[lev][0]->nGrow());
430 cmesh_exdata.setVal(0.0, 0, 1, mesh_data[lev][0]->nGrow());
431 cmesh_exdata.copy(*mesh_data[lev-1][0], 0, 0, 1, 1, 1);
432 cmesh_exdata.FillBoundary(gm.periodicity());
435 mesh_data[lev][1]->nComp(),
436 mesh_data[lev][1]->nGrow());
437 cmesh_eydata.setVal(0.0, 0, 1, mesh_data[lev][1]->nGrow());
438 cmesh_eydata.copy(*mesh_data[lev-1][1], 0, 0, 1, 1, 1);
439 cmesh_eydata.FillBoundary(gm.periodicity());
442 mesh_data[lev][2]->nComp(),
443 mesh_data[lev][2]->nGrow());
444 cmesh_ezdata.setVal(0.0, 0, 1, mesh_data[lev][2]->nGrow());
445 cmesh_ezdata.copy(*mesh_data[lev-1][2], 0, 0, 1, 1, 1);
446 cmesh_ezdata.FillBoundary(gm.periodicity());
448 for (
size_t ip = lBegin; ip < lEnd; ++ip) {
450 const int grid = this->Grid[ip];
460 const typename PLayout::basefab_t& mfab = (*layout_p->getLevelMask(lev))[grid];
464 for (
int ii = 0; ii < 3; ++ii) {
465 lxyz[ii] = ( this->
R[ip](ii) - plo[ii] ) * inv_fdx[ii] + 0.5;
467 wxyz_hi[ii] = lxyz[ii] - ijk[ii];
468 wxyz_lo[ii] = 1.0 - wxyz_hi[ii];
475 bool use_coarse =
false;
482 for (
int ii = 0; ii < 3; ++ii) {
483 lxyz[ii] = ( this->
R[ip](ii) - cplo[ii] ) * inv_cdx[ii] + 0.5;
485 wxyz_hi[ii] = lxyz[ii] - ijk[ii];
486 wxyz_lo[ii] = 1.0 - wxyz_hi[ii];
501 pa[ip](0) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * cexfab(i1) +
502 wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * cexfab(i3) +
503 wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * cexfab(i5) +
504 wxyz_hi[0] * wxyz_hi[1] * wxyz_lo[2] * cexfab(i7) +
505 wxyz_lo[0] * wxyz_lo[1] * wxyz_hi[2] * cexfab(i2) +
506 wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * cexfab(i4) +
507 wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * cexfab(i6) +
508 wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * cexfab(i8);
510 pa[ip](1) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * ceyfab(i1) +
511 wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * ceyfab(i3) +
512 wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * ceyfab(i5) +
513 wxyz_hi[0] * wxyz_hi[1] * wxyz_lo[2] * ceyfab(i7) +
514 wxyz_lo[0] * wxyz_lo[1] * wxyz_hi[2] * ceyfab(i2) +
515 wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * ceyfab(i4) +
516 wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * ceyfab(i6) +
517 wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * ceyfab(i8);
519 pa[ip](2) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * cezfab(i1) +
520 wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * cezfab(i3) +
521 wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * cezfab(i5) +
522 wxyz_hi[0] * wxyz_hi[1] * wxyz_lo[2] * cezfab(i7) +
523 wxyz_lo[0] * wxyz_lo[1] * wxyz_hi[2] * cezfab(i2) +
524 wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * cezfab(i4) +
525 wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * cezfab(i6) +
526 wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * cezfab(i8);
528 pa[ip](0) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * exfab(i1) +
529 wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * exfab(i3) +
530 wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * exfab(i5) +
531 wxyz_hi[0] * wxyz_hi[1] * wxyz_lo[2] * exfab(i7) +
532 wxyz_lo[0] * wxyz_lo[1] * wxyz_hi[2] * exfab(i2) +
533 wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * exfab(i4) +
534 wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * exfab(i6) +
535 wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * exfab(i8);
537 pa[ip](1) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * eyfab(i1) +
538 wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * eyfab(i3) +
539 wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * eyfab(i5) +
540 wxyz_hi[0] * wxyz_hi[1] * wxyz_lo[2] * eyfab(i7) +
541 wxyz_lo[0] * wxyz_lo[1] * wxyz_hi[2] * eyfab(i2) +
542 wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * eyfab(i4) +
543 wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * eyfab(i6) +
544 wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * eyfab(i8);
546 pa[ip](2) = wxyz_lo[0] * wxyz_lo[1] * wxyz_lo[2] * ezfab(i1) +
547 wxyz_lo[0] * wxyz_hi[1] * wxyz_lo[2] * ezfab(i3) +
548 wxyz_hi[0] * wxyz_lo[1] * wxyz_lo[2] * ezfab(i5) +
549 wxyz_hi[0] * wxyz_hi[1] * wxyz_lo[2] * ezfab(i7) +
550 wxyz_lo[0] * wxyz_lo[1] * wxyz_hi[2] * ezfab(i2) +
551 wxyz_lo[0] * wxyz_hi[1] * wxyz_hi[2] * ezfab(i4) +
552 wxyz_hi[0] * wxyz_lo[1] * wxyz_hi[2] * ezfab(i6) +
553 wxyz_hi[0] * wxyz_hi[1] * wxyz_hi[2] * ezfab(i8);
void InterpolateMultiLevelFort(ParticleAttrib< AType > &pa, AmrVectorFieldContainer_t &mesh_data, int lev)
PLayout::AmrGrid_t AmrGrid_t
void AssignDensityFort(ParticleAttrib< AType > &pa, AmrScalarFieldContainer_t &mf_to_be_filled, int lev_min, int ncomp, int finest_level, const ParticleAttrib< int > &pbin, int bin=-1) const
void InterpolateSingleLevelFort(ParticleAttrib< AType > &pa, AmrVectorField_t &mesh_data, int lev)
AmrParticleBase< PLayout >::AmrField_t AmrField_t
The base class for all OPAL exceptions.
AmrParticleBase< PLayout >::AmrScalarFieldContainer_t AmrScalarFieldContainer_t
amrex::IntVect AmrIntVect_t
AmrParticleBase< PLayout >::AmrVectorFieldContainer_t AmrVectorFieldContainer_t
bool scatter(Communicate &, InputIterator, InputIterator, RandomIterator, int *, int *, const ScatterOp &)
void AssignCellDensitySingleLevelFort(ParticleAttrib< AType > &pa, AmrField_t &mf, int level, const ParticleAttrib< int > &pbin, int bin=-1, int ncomp=1, int particle_lvl_offset=0) const
PLayout::AmrGeometry_t AmrGeometry_t
IpplTimings::TimerRef AssignDensityTimer_m
PLayout::AmrReal_t AmrReal_t
static void startTimer(TimerRef t)
void scatter(ParticleAttrib< FT > &attrib, AmrScalarFieldContainer_t &f, ParticleAttrib< Vektor< PT, Dim > > &pp, int lbase, int lfine, const ParticleAttrib< int > &pbin, int bin=-1)
void InterpolateFort(ParticleAttrib< AType > &pa, AmrVectorFieldContainer_t &mesh_data, int lev_min, int lev_max)
PLayout::AmrProcMap_t AmrProcMap_t
AmrParticleBase< PLayout >::AmrVectorField_t AmrVectorField_t
amrex::FArrayBox FArrayBox_t
static TimerRef getTimer(const char *nm)
static void stopTimer(TimerRef t)
PLayout::AmrIntVect_t AmrIntVect_t
void gather(ParticleAttrib< FT > &attrib, AmrVectorFieldContainer_t &f, ParticleAttrib< Vektor< PT, Dim > > &pp, int lbase, int lfine)
amr::AmrField_t AmrField_t