1 #ifndef BoxLibLayout_HPP
2 #define BoxLibLayout_HPP
12 template <
class T,
unsigned Dim>
16 template <
class T,
unsigned Dim>
20 template<
class T,
unsigned Dim>
38 int nGridPoints =
std::ceil( std::cbrt( nProcs ) ) * maxGridSize;
44 template<
class T,
unsigned Dim>
47 ParGDB(layout_p->m_geom,
59 template<
class T,
unsigned Dim>
69 template<
class T,
unsigned Dim>
74 ParGDB(geom, dmap, ba),
79 template<
class T,
unsigned Dim>
85 ParGDB(geom, dmap, ba, rr),
90 template<
class T,
unsigned Dim>
94 int nGridPoints = this->m_geom[0].Domain().length(0);
95 int maxGridSize = this->m_ba[0][0].length(0);
97 this->initBaseBox_m(nGridPoints, maxGridSize, dh);
101 template <
class T,
unsigned Dim>
104 static bool isCalled =
false;
112 if ( ratio.size() !=
Dim ) {
114 "Length " + std::to_string(ratio.size()) +
115 " != " + std::to_string(
Dim));
118 for (
unsigned int i = 0; i <
Dim; ++i) {
119 if ( ratio[i] <= 0.0 ) {
121 "The ratio has to be larger than zero.");
124 lowerBound[i] *= ratio[i];
125 upperBound[i] *= ratio[i];
130 template<
class T,
unsigned Dim>
137 throw OpalException(
"BoxLibLayout::update(IpplParticleBase, ParticleAttrib) ",
138 "Wrong update method called.");
144 template<
class T,
unsigned Dim>
146 int lev_min,
int lev_max,
bool isRegrid)
149 if ( !PData.isForbidTransform() ) {
153 PData.domainMapping();
161 int theEffectiveFinestLevel = this->finestLevel();
162 while (!this->LevelDefined(theEffectiveFinestLevel)) {
163 theEffectiveFinestLevel--;
167 lev_max = theEffectiveFinestLevel;
168 else if ( lev_max > theEffectiveFinestLevel )
169 lev_max = theEffectiveFinestLevel;
172 size_t LocalNum = PData.getLocalNum();
174 auto& LocalNumPerLevel = PData.getLocalNumPerLevel();
176 if ( LocalNum != LocalNumPerLevel.getLocalNumAllLevel() )
178 "Local #particles disagrees with sum over levels");
180 std::multimap<unsigned, unsigned> p2n;
182 int *msgsend =
new int[N];
183 std::fill(msgsend, msgsend+N, 0);
184 int *msgrecv =
new int[N];
185 std::fill(msgrecv, msgrecv+N, 0);
188 size_t lBegin = LocalNumPerLevel.begin(lev_min);
189 size_t lEnd = LocalNumPerLevel.end(lev_max);
202 for (
unsigned int ip = lBegin; ip < lEnd; ++ip) {
204 const size_t& lold = PData.Level[ip];
213 locateParticle(PData, ip, lev_min, lev_max, nGrow);
217 const size_t& lnew = PData.Level[ip];
219 const unsigned int who = ParticleDistributionMap(lnew)[PData.Grid[ip]];
221 --LocalNumPerLevel[lold];
226 p2n.insert(std::pair<unsigned, unsigned>(who, ip));
232 ++LocalNumPerLevel[lnew];
237 MPI_Allreduce(msgsend, msgrecv, N, MPI_INT, MPI_SUM,
Ippl::getComm());
243 Format *format = PData.getFormat();
245 std::vector<MPI_Request> requests;
246 std::vector<MsgBuffer*> buffers;
251 unsigned cur_destination = i->first;
255 for (; i!=p2n.end() && i->first == cur_destination; ++i)
258 PData.putMessage(msg, i->second);
259 PData.destroy(1, i->second);
265 cur_destination, tag);
268 requests.push_back(request);
269 buffers.push_back(msgbuf);
274 if ( LocalNum < PData.getDestroyNum() )
276 "Rank " + std::to_string(myN) +
277 " can't destroy more particles than possessed.");
279 LocalNum -= PData.getDestroyNum();
280 PData.performDestroy();
283 for (
int lev = lev_min; lev <= lev_max; ++lev) {
284 if ( LocalNumPerLevel[lev] < 0 )
286 "Negative particle level count.");
290 for (
int k = 0; k<msgrecv[myN]; ++k)
295 MsgBuffer recvbuf(format, buffer, bufsize);
303 size_t pBeginIdx = LocalNum;
305 LocalNum += PData.getSingleMessage(*msg);
307 size_t pEndIdx = LocalNum;
309 for (
size_t idx = pBeginIdx; idx < pEndIdx; ++idx)
310 ++LocalNumPerLevel[ PData.Level[idx] ];
318 MPI_Waitall(requests.size(), &(requests[0]), MPI_STATUSES_IGNORE);
319 for (
unsigned int j = 0; j<buffers.size(); ++j)
334 MPI_Allreduce(&LocalNum, &TotalNum, 1, MPI_INT, MPI_SUM,
Ippl::getComm());
337 PData.setTotalNum(TotalNum);
338 PData.setLocalNum(LocalNum);
341 if ( LocalNum != LocalNumPerLevel.getLocalNumAllLevel() )
343 "Local #particles disagrees with sum over levels");
345 if ( !PData.isForbidTransform() ) {
347 PData.domainMapping(
true);
354 template <
class T,
unsigned Dim>
357 const unsigned int ip,
360 return Index(p.R[ip], lev);
364 template <
class T,
unsigned Dim>
372 D_TERM(iv[0]=
floor((R[0]-geom.ProbLo(0))/geom.CellSize(0));,
373 iv[1]=
floor((R[1]-geom.ProbLo(1))/geom.CellSize(1));,
374 iv[2]=
floor((R[2]-geom.ProbLo(2))/geom.CellSize(2)););
376 iv += geom.Domain().smallEnd();
382 template <
class T,
unsigned Dim>
389 if ( lev >= (
int)masks_m.size() )
390 masks_m.resize(lev + 1);
392 masks_m[lev].reset(
new mask_t(ParticleBoxArray(lev),
393 ParticleDistributionMap(lev), 1, 1));
395 masks_m[lev]->setVal(1, 1);
397 mask_t tmp_mask(ParticleBoxArray(lev),
398 ParticleDistributionMap(lev),
401 tmp_mask.setVal(0, ncells);
403 tmp_mask.BuildMask(Geom(lev).Domain(), Geom(lev).periodicity(),
404 covered, notcovered, physbnd, interior);
406 tmp_mask.FillBoundary(Geom(lev).periodicity());
408 for (amrex::MFIter mfi(tmp_mask); mfi.isValid(); ++mfi) {
409 const AmrBox_t& bx = mfi.validbox();
410 const int* lo = bx.loVect();
411 const int* hi = bx.hiVect();
416 for (
int i = lo[0]; i <= hi[0]; ++i) {
417 for (
int j = lo[1]; j <= hi[1]; ++j) {
418 for (
int k = lo[2]; k <= hi[2]; ++k) {
421 for (
int ii = i - ncells; ii <= i + ncells; ++ii) {
422 for (
int jj = j - ncells; jj <= j + ncells; ++jj) {
423 for (
int kk = k - ncells; kk <= k + ncells; ++kk) {
439 masks_m[lev]->FillBoundary(Geom(lev).periodicity());
443 template <
class T,
unsigned Dim>
445 assert(lev < (
int)masks_m.size());
446 masks_m[lev].reset(
nullptr);
450 template <
class T,
unsigned Dim>
451 const std::unique_ptr<typename BoxLibLayout<T, Dim>::mask_t>&
453 if ( lev >= (
int)masks_m.size() ) {
455 "Unable to access level " + std::to_string(lev) +
".");
509 template <
class T,
unsigned Dim>
511 const unsigned int ip,
518 lev_max = finestLevel();
520 BL_ASSERT(lev_max <= finestLevel());
522 BL_ASSERT(nGrow == 0 || (nGrow >= 0 && lev_min == lev_max));
524 std::vector< std::pair<int, AmrBox_t> > isects;
526 for (
int lev = lev_max; lev >= lev_min; lev--)
529 const AmrGrid_t& ba = ParticleBoxArray(lev);
530 BL_ASSERT(ba.ixType().cellCentered());
532 if (lev == (
int)p.Level[ip]) {
534 if (0 <= p.Grid[ip] && p.Grid[ip] < ba.size())
536 const AmrBox_t& bx = ba.getCellCenteredBox(p.Grid[ip]);
537 const AmrBox_t& gbx = amrex::grow(bx,nGrow);
538 if (gbx.contains(iv))
549 ba.intersections(
AmrBox_t(iv, iv), isects,
true, nGrow);
554 p.Grid[ip] = isects[0].first;
566 template <
class T,
unsigned Dim>
568 const unsigned int ip,
572 if (!Geom(0).isAnyPeriodic())
return false;
575 lev_max = finestLevel();
577 BL_ASSERT(lev_max <= finestLevel());
583 if (PeriodicShift(R))
585 std::vector< std::pair<int, AmrBox_t> > isects;
587 for (
int lev = lev_max; lev >= lev_min; lev--)
590 const AmrGrid_t& ba = ParticleBoxArray(lev);
592 ba.intersections(
AmrBox_t(iv,iv),isects,
true,0);
596 D_TERM(p.R[ip][0] = R[0];,
601 p.Grid[ip] = isects[0].first;
614 template <
class T,
unsigned Dim>
624 const AmrBox_t& dmn = geom.Domain();
626 bool shifted =
false;
628 for (
int i = 0; i < AMREX_SPACEDIM; i++)
630 if (!geom.isPeriodic(i))
continue;
632 if (iv[i] > dmn.bigEnd(i))
634 if (R[i] == geom.ProbHi(i))
640 R[i] += .125*geom.CellSize(i);
642 R[i] -= geom.ProbLength(i);
644 if (R[i] <= geom.ProbLo(i))
648 R[i] += .125*geom.CellSize(i);
650 BL_ASSERT(R[i] >= geom.ProbLo(i));
654 else if (iv[i] < dmn.smallEnd(i))
656 if (R[i] == geom.ProbLo(i))
662 R[i] -= .125*geom.CellSize(i);
664 R[i] += geom.ProbLength(i);
666 if (R[i] >= geom.ProbHi(i))
670 R[i] -= .125*geom.CellSize(i);
672 BL_ASSERT(R[i] <= geom.ProbHi(i));
685 template <
class T,
unsigned Dim>
688 const unsigned int ip,
689 int lev_min,
int lev_max,
int nGrow)
const
691 bool outside = D_TERM( p.R[ip](0) < AmrGeometry_t::ProbLo(0)
692 || p.R[ip](0) >= AmrGeometry_t::ProbHi(0),
693 || p.R[ip](1) < AmrGeometry_t::ProbLo(1)
694 || p.R[ip](1) >= AmrGeometry_t::ProbHi(1),
695 || p.R[ip](2) < AmrGeometry_t::ProbLo(2)
696 || p.R[ip](2) >= AmrGeometry_t::ProbHi(2));
698 bool success =
false;
703 success = EnforcePeriodicWhere(p, ip, lev_min, lev_max);
704 if (!success && lev_min == 0)
714 "We're losing particles although we shouldn't");
720 success = Where(p, ip, lev_min, lev_max);
725 success = (nGrow > 0) && Where(p, ip, lev_min, lev_min, nGrow);
730 std::stringstream ss;
731 ss <<
"Invalid particle with ID " << ip <<
" at position " << p.R[ip] <<
".";
732 throw OpalException(
"BoxLibLayout::locateParticle()", ss.str());
738 template <
class T,
unsigned Dim>
740 return level <= this->maxLevel_m && !m_ba[level].empty() && !m_dmap[level].empty();
744 template <
class T,
unsigned Dim>
746 return this->finestLevel_m;
750 template <
class T,
unsigned Dim>
752 return this->maxLevel_m;
756 template <
class T,
unsigned Dim>
760 return refRatio_m[level];
764 template <
class T,
unsigned Dim>
767 for (
int n = 0;
n<AMREX_SPACEDIM;
n++)
768 maxval =
std::max(maxval, refRatio_m[level][
n]);
773 template <
class T,
unsigned Dim>
780 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
782 assert(lowerBound[d] < 0);
783 assert(upperBound[d] > 0);
785 real_box.setLo(d, lowerBound[d] * (1.0 + dh));
786 real_box.setHi(d, upperBound[d] * (1.0 + dh));
789 AmrGeometry_t::ProbDomain(real_box);
793 AmrIntVect_t domain_hi(nGridPoints - 1, nGridPoints - 1, nGridPoints - 1);
794 const AmrBox_t domain(domain_lo, domain_hi);
800 int is_per[AMREX_SPACEDIM];
801 for (
int i = 0; i < AMREX_SPACEDIM; i++)
805 geom.define(domain, &real_box, coord, is_per);
810 ba.maxSize(maxGridSize);
816 this->m_geom.resize(1);
817 this->m_geom[0] = geom;
819 this->m_dmap.resize(1);
820 this->m_dmap[0] = dmap;
822 this->m_ba.resize(1);
825 this->m_nlevels = ba.size();
void update(IpplParticleBase< BoxLibLayout< T, Dim > > &PData, const ParticleAttrib< char > *canSwap=0)
AmrIntVect_t refRatio(int level) const
const std::unique_ptr< mask_t > & getLevelMask(int lev) const
void locateParticle(AmrParticleBase< BoxLibLayout< T, Dim > > &p, const unsigned int ip, int lev_min, int lev_max, int nGrow) const
amrex::BaseFab< int > basefab_t
amr::AmrIntVect_t AmrIntVect_t
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
The base class for all OPAL exceptions.
amr::AmrGeomContainer_t AmrGeomContainer_t
PETE_TUTree< FnCeil, typename T::PETE_Expr_t > ceil(const PETE_Expr< T > &l)
virtual MPI_Request raw_isend(void *, int, int, int)
virtual int raw_probe_receive(char *&, int &, int &)
amr::AmrProcMapContainer_t AmrProcMapContainer_t
int next_tag(int t, int s=1000)
amrex::FabArray< basefab_t > mask_t
amr::AmrDomain_t AmrDomain_t
amr::AmrGridContainer_t AmrGridContainer_t
void setDomainRatio(const std::vector< double > &ratio)
bool Where(AmrParticleBase< BoxLibLayout< T, Dim > > &p, const unsigned int ip, int lev_min=0, int lev_max=-1, int nGrow=0) const
bool EnforcePeriodicWhere(AmrParticleBase< BoxLibLayout< T, Dim > > &prt, const unsigned int ip, int lev_min=0, int lev_max=-1) const
int MaxRefRatio(int level) const
void clearLevelMask(int lev)
AmrIntVect_t Index(AmrParticleBase< BoxLibLayout< T, Dim > > &p, const unsigned int ip, int level) const
Vektor< double, 3 > Vector_t
amr::AmrGeometry_t AmrGeometry_t
static MPI_Comm getComm()
void buildLevelMask(int lev, const int ncells=1)
int maxLevel_m
Maximum level allowed.
bool LevelDefined(int level) const
amr::AmrProcMap_t AmrProcMap_t
bool PeriodicShift(SingleParticlePos_t R) const
amr::AmrIntArray_t AmrIntArray_t
void initBaseBox_m(int nGridPoints, int maxGridSize, double dh=0.04)
void setBoundingBox(double dh)
std::string::iterator iterator
AmrIntVectContainer_t refRatio_m
static Communicate * Comm
#define P_SPATIAL_TRANSFER_TAG
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)