BoxLibLayout.hpp 25.4 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6
#ifndef BoxLibLayout_HPP
#define BoxLibLayout_HPP

#include "BoxLibLayout.h"

#include "Message/Formatter.h"
7
#include "Utilities/OpalException.h"
frey_m's avatar
frey_m committed
8

9
#include <cassert>
frey_m's avatar
frey_m committed
10 11
#include <cmath>

frey_m's avatar
frey_m committed
12
template <class T, unsigned Dim>
13
Vector_t BoxLibLayout<T, Dim>::lowerBound = - Vector_t(1.0, 1.0, 1.0);
frey_m's avatar
frey_m committed
14 15 16


template <class T, unsigned Dim>
17
Vector_t BoxLibLayout<T, Dim>::upperBound = Vector_t(1.0, 1.0, 1.0);
frey_m's avatar
frey_m committed
18 19


20
template<class T, unsigned Dim>
frey_m's avatar
frey_m committed
21
BoxLibLayout<T, Dim>::BoxLibLayout()
frey_m's avatar
frey_m committed
22 23 24
    : ParticleAmrLayout<T, Dim>(),
      ParGDB(),
      refRatio_m(0)
frey_m's avatar
frey_m committed
25
{
frey_m's avatar
frey_m committed
26
    /* FIXME There might be a better solution
kraus's avatar
kraus committed
27 28
     *
     *
frey_m's avatar
frey_m committed
29 30
     * Figure out the number of grid points in each direction
     * such that all processes have some data at the beginning
kraus's avatar
kraus committed
31
     *
frey_m's avatar
frey_m committed
32
     * ( nGridPoints / maxGridSize ) ^3 = max. #procs
kraus's avatar
kraus committed
33
     *
frey_m's avatar
frey_m committed
34 35 36
     */
    int nProcs = Ippl::getNodes();
    int maxGridSize = 16;
kraus's avatar
kraus committed
37

frey_m's avatar
frey_m committed
38
    int nGridPoints = std::ceil( std::cbrt( nProcs ) ) * maxGridSize;
kraus's avatar
kraus committed
39

frey_m's avatar
frey_m committed
40
    this->initBaseBox_m(nGridPoints, maxGridSize);
frey_m's avatar
frey_m committed
41
}
frey_m's avatar
frey_m committed
42 43


frey_m's avatar
frey_m committed
44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
template<class T, unsigned Dim>
BoxLibLayout<T, Dim>::BoxLibLayout(const BoxLibLayout* layout_p)
    : ParticleAmrLayout<T, Dim>(),
      ParGDB(layout_p->m_geom,
             layout_p->m_dmap,
             layout_p->m_ba,
             layout_p->m_rr)
{
    this->maxLevel_m = layout_p->maxLevel_m;
    refRatio_m.resize(layout_p->m_geom.size()-1);
    for (int i = 0; i < refRatio_m.size(); ++i)
        refRatio_m[i] = layout_p->refRatio(i);
}


59
template<class T, unsigned Dim>
frey_m's avatar
frey_m committed
60
BoxLibLayout<T, Dim>::BoxLibLayout(int nGridPoints, int maxGridSize)
frey_m's avatar
frey_m committed
61 62 63
    : ParticleAmrLayout<T, Dim>(),
      ParGDB(),
      refRatio_m(0)
64
{
frey_m's avatar
frey_m committed
65
    this->initBaseBox_m(nGridPoints, maxGridSize);
66 67 68
}


frey_m's avatar
frey_m committed
69
template<class T, unsigned Dim>
frey_m's avatar
frey_m committed
70 71 72
BoxLibLayout<T, Dim>::BoxLibLayout(const AmrGeometry_t &geom,
                                   const AmrProcMap_t &dmap,
                                   const AmrGrid_t &ba)
frey_m's avatar
frey_m committed
73 74 75
    : ParticleAmrLayout<T, Dim>(),
      ParGDB(geom, dmap, ba),
      refRatio_m(0)
frey_m's avatar
frey_m committed
76 77 78 79
{ }


template<class T, unsigned Dim>
frey_m's avatar
frey_m committed
80 81 82 83
BoxLibLayout<T, Dim>::BoxLibLayout(const AmrGeomContainer_t &geom,
                                   const AmrProcMapContainer_t &dmap,
                                   const AmrGridContainer_t &ba,
                                   const AmrIntArray_t &rr)
frey_m's avatar
frey_m committed
84 85 86
    : ParticleAmrLayout<T, Dim>(),
      ParGDB(geom, dmap, ba, rr),
      refRatio_m(0)
frey_m's avatar
frey_m committed
87 88 89
{ }


frey_m's avatar
frey_m committed
90 91
template<class T, unsigned Dim>
void BoxLibLayout<T, Dim>::setBoundingBox(double dh) {
kraus's avatar
kraus committed
92

frey_m's avatar
frey_m committed
93 94 95
    // cubic box
    int nGridPoints = this->m_geom[0].Domain().length(0);
    int maxGridSize = this->m_ba[0][0].length(0);
kraus's avatar
kraus committed
96

frey_m's avatar
frey_m committed
97 98 99 100
    this->initBaseBox_m(nGridPoints, maxGridSize, dh);
}


101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129
template <class T, unsigned Dim>
void BoxLibLayout<T, Dim>::setDomainRatio(const std::vector<double>& ratio) {

    static bool isCalled = false;

    if ( isCalled ) {
        return;
    }

    isCalled = true;

    if ( ratio.size() != Dim ) {
        throw OpalException("BoxLibLayout::setDomainRatio() ",
                            "Length " + std::to_string(ratio.size()) +
                            " != " + std::to_string(Dim));
    }

    for (unsigned int i = 0; i < Dim; ++i) {
        if ( ratio[i] <= 0.0 ) {
            throw OpalException("BoxLibLayout::setDomainRatio() ",
                                "The ratio has to be larger than zero.");
        }

        lowerBound[i] *= ratio[i];
        upperBound[i] *= ratio[i];
    }
}


frey_m's avatar
frey_m committed
130 131 132
template<class T, unsigned Dim>
void BoxLibLayout<T, Dim>::update(IpplParticleBase< BoxLibLayout<T,Dim> >& PData,
                                  const ParticleAttrib<char>* canSwap)
133
{
134
    /* Exit since we need AmrParticleBase with grids and levels for particles for this layout
frey_m's avatar
frey_m committed
135 136
     * if IpplParticleBase is used something went wrong
     */
137 138
    throw OpalException("BoxLibLayout::update(IpplParticleBase, ParticleAttrib) ",
                        "Wrong update method called.");
139
}
frey_m's avatar
frey_m committed
140

frey_m's avatar
frey_m committed
141

frey_m's avatar
frey_m committed
142
// // Function from AMReX adjusted to work with Ippl AmrParticleBase class
frey_m's avatar
frey_m committed
143 144
// // redistribute the particles using BoxLibs ParGDB class to determine where particle should go
template<class T, unsigned Dim>
145
void BoxLibLayout<T, Dim>::update(AmrParticleBase< BoxLibLayout<T,Dim> >& PData,
146
                                  int lev_min, int lev_max, bool isRegrid)
frey_m's avatar
frey_m committed
147
{
148
    // in order to avoid transforms when already done
frey_m's avatar
frey_m committed
149
    if ( !PData.isForbidTransform() ) {
150 151 152
        /* we need to update on Amr domain + boosted frame (Lorentz transform,
         * has to be undone at end of function
         */
frey_m's avatar
frey_m committed
153
        PData.domainMapping();
154
    }
kraus's avatar
kraus committed
155

frey_m's avatar
frey_m committed
156 157 158 159
    int nGrow = 0;

    unsigned N = Ippl::getNodes();
    unsigned myN = Ippl::myNode();
kraus's avatar
kraus committed
160

frey_m's avatar
frey_m committed
161
    int theEffectiveFinestLevel = this->finestLevel();
frey_m's avatar
frey_m committed
162
    while (!this->LevelDefined(theEffectiveFinestLevel)) {
frey_m's avatar
frey_m committed
163
        theEffectiveFinestLevel--;
frey_m's avatar
frey_m committed
164
    }
kraus's avatar
kraus committed
165

frey_m's avatar
frey_m committed
166 167
    if (lev_max == -1)
        lev_max = theEffectiveFinestLevel;
frey_m's avatar
frey_m committed
168 169
    else if ( lev_max > theEffectiveFinestLevel )
        lev_max = theEffectiveFinestLevel;
kraus's avatar
kraus committed
170

171
    //loop trough the particles and assign the grid and level where each particle belongs
frey_m's avatar
frey_m committed
172
    size_t LocalNum = PData.getLocalNum();
kraus's avatar
kraus committed
173

frey_m's avatar
frey_m committed
174
    auto& LocalNumPerLevel = PData.getLocalNumPerLevel();
kraus's avatar
kraus committed
175

frey_m's avatar
frey_m committed
176 177 178
    if ( LocalNum != LocalNumPerLevel.getLocalNumAllLevel() )
        throw OpalException("BoxLibLayout::update()",
                            "Local #particles disagrees with sum over levels");
kraus's avatar
kraus committed
179 180

    std::multimap<unsigned, unsigned> p2n; //node ID, particle
frey_m's avatar
frey_m committed
181 182 183 184 185 186 187

    int *msgsend = new int[N];
    std::fill(msgsend, msgsend+N, 0);
    int *msgrecv = new int[N];
    std::fill(msgrecv, msgrecv+N, 0);

    unsigned sent = 0;
frey_m's avatar
frey_m committed
188 189
    size_t lBegin = LocalNumPerLevel.begin(lev_min);
    size_t lEnd   = LocalNumPerLevel.end(lev_max);
kraus's avatar
kraus committed
190

191 192 193 194 195 196 197 198 199
    /* in the case of regrid we might lose a level
     * and therefore the level counter is invalidated
     */
    if ( isRegrid ) {
        lBegin = 0;
        lEnd   = LocalNum;
    }


frey_m's avatar
frey_m committed
200 201
    //loop trough particles and assign grid and level to each particle
    //if particle doesn't belong to this process save the index of the particle to be sent
frey_m's avatar
frey_m committed
202
    for (unsigned int ip = lBegin; ip < lEnd; ++ip) {
frey_m's avatar
frey_m committed
203 204
        // old level
        const size_t& lold = PData.Level[ip];
kraus's avatar
kraus committed
205

frey_m's avatar
frey_m committed
206
//         /*
frey_m's avatar
frey_m committed
207
//          * AMReX sets m_grid = -1 and m_lev = -1
frey_m's avatar
frey_m committed
208 209 210
//          */
//         PData.Level[ip] = -1;
//         PData.Grid[ip] = -1;
kraus's avatar
kraus committed
211

frey_m's avatar
frey_m committed
212
        //check to which level and grid the particle belongs to
frey_m's avatar
frey_m committed
213
        locateParticle(PData, ip, lev_min, lev_max, nGrow);
kraus's avatar
kraus committed
214

frey_m's avatar
frey_m committed
215 216 217
        // The owner of the particle is the CPU owning the finest grid
        // in state data that contains the particle.
        const size_t& lnew = PData.Level[ip];
kraus's avatar
kraus committed
218

frey_m's avatar
frey_m committed
219
        const unsigned int who = ParticleDistributionMap(lnew)[PData.Grid[ip]];
kraus's avatar
kraus committed
220

frey_m's avatar
frey_m committed
221
        --LocalNumPerLevel[lold];
kraus's avatar
kraus committed
222

frey_m's avatar
frey_m committed
223 224 225 226 227
        if (who != myN) {
            // we lost the particle to another process
            msgsend[who] = 1;
            p2n.insert(std::pair<unsigned, unsigned>(who, ip));
            sent++;
228
        } else {
frey_m's avatar
frey_m committed
229 230
            /* if we still own the particle it may have moved to
             * another level
231
             */
frey_m's avatar
frey_m committed
232
            ++LocalNumPerLevel[lnew];
frey_m's avatar
frey_m committed
233
        }
frey_m's avatar
frey_m committed
234
    }
frey_m's avatar
frey_m committed
235 236

    //reduce message count so every node knows how many messages to receive
237
    allreduce(msgsend.data(), msgrecv.data(), N, std::plus<int>());
frey_m's avatar
frey_m committed
238 239 240 241 242 243 244 245 246 247 248 249 250 251

    int tag = Ippl::Comm->next_tag(P_SPATIAL_TRANSFER_TAG,P_LAYOUT_CYCLE);

    typename std::multimap<unsigned, unsigned>::iterator i = p2n.begin();

    Format *format = PData.getFormat();

    std::vector<MPI_Request> requests;
    std::vector<MsgBuffer*> buffers;

    //create a message and send particles to nodes they belong to
    while (i!=p2n.end())
    {
        unsigned cur_destination = i->first;
kraus's avatar
kraus committed
252

frey_m's avatar
frey_m committed
253 254 255 256 257 258 259 260 261
        MsgBuffer *msgbuf = new MsgBuffer(format, p2n.count(i->first));

        for (; i!=p2n.end() && i->first == cur_destination; ++i)
        {
            Message msg;
            PData.putMessage(msg, i->second);
            PData.destroy(1, i->second);
            msgbuf->add(&msg);
        }
kraus's avatar
kraus committed
262 263 264

        MPI_Request request = Ippl::Comm->raw_isend(msgbuf->getBuffer(),
                                                    msgbuf->getSize(),
frey_m's avatar
frey_m committed
265 266 267 268 269
                                                    cur_destination, tag);

        //remember request and buffer so we can delete them later
        requests.push_back(request);
        buffers.push_back(msgbuf);
kraus's avatar
kraus committed
270

frey_m's avatar
frey_m committed
271
    }
kraus's avatar
kraus committed
272

frey_m's avatar
frey_m committed
273
    //destroy the particles that are sent to other domains
frey_m's avatar
frey_m committed
274
    if ( LocalNum < PData.getDestroyNum() )
frey_m's avatar
frey_m committed
275 276 277
        throw OpalException("BoxLibLayout::update()",
                            "Rank " + std::to_string(myN) +
                            " can't destroy more particles than possessed.");
frey_m's avatar
frey_m committed
278 279 280 281
    else {
        LocalNum -= PData.getDestroyNum();  // update local num
        PData.performDestroy();
    }
kraus's avatar
kraus committed
282

frey_m's avatar
frey_m committed
283 284 285 286 287
    for (int lev = lev_min; lev <= lev_max; ++lev) {
        if ( LocalNumPerLevel[lev] < 0 )
            throw OpalException("BoxLibLayout::update()",
                                "Negative particle level count.");
    }
kraus's avatar
kraus committed
288

frey_m's avatar
frey_m committed
289 290 291 292 293 294 295
    //receive new particles
    for (int k = 0; k<msgrecv[myN]; ++k)
    {
        int node = Communicate::COMM_ANY_NODE;
        char *buffer = 0;
        int bufsize = Ippl::Comm->raw_probe_receive(buffer, node, tag);
        MsgBuffer recvbuf(format, buffer, bufsize);
kraus's avatar
kraus committed
296

frey_m's avatar
frey_m committed
297 298 299
        Message *msg = recvbuf.get();
        while (msg != 0)
            {
frey_m's avatar
frey_m committed
300 301 302 303
                /* pBeginIdx is the start index of the new particle data
                 * pEndIdx is the last index of the new particle data
                 */
                size_t pBeginIdx = LocalNum;
kraus's avatar
kraus committed
304

frey_m's avatar
frey_m committed
305
                LocalNum += PData.getSingleMessage(*msg);
kraus's avatar
kraus committed
306

frey_m's avatar
frey_m committed
307
                size_t pEndIdx = LocalNum;
kraus's avatar
kraus committed
308

frey_m's avatar
frey_m committed
309 310
                for (size_t idx = pBeginIdx; idx < pEndIdx; ++idx)
                    ++LocalNumPerLevel[ PData.Level[idx] ];
kraus's avatar
kraus committed
311

frey_m's avatar
frey_m committed
312 313
                delete msg;
                msg = recvbuf.get();
kraus's avatar
kraus committed
314
            }
frey_m's avatar
frey_m committed
315
    }
kraus's avatar
kraus committed
316

frey_m's avatar
frey_m committed
317 318 319 320 321 322
    //wait for communication to finish and clean up buffers
    MPI_Waitall(requests.size(), &(requests[0]), MPI_STATUSES_IGNORE);
    for (unsigned int j = 0; j<buffers.size(); ++j)
    {
        delete buffers[j];
    }
kraus's avatar
kraus committed
323

frey_m's avatar
frey_m committed
324 325 326 327 328 329 330 331 332 333
    delete[] msgsend;
    delete[] msgrecv;
    delete format;

    // there is extra work to do if there are multipple nodes, to distribute
    // the particle layout data to all nodes
    //TODO: do we need info on how many particles are on each node?

    //save how many total particles we have
    size_t TotalNum = 0;
334
    allreduce(&LocalNum, &TotalNum, 1, std::plus<size_t>());
frey_m's avatar
frey_m committed
335 336

    // update our particle number counts
frey_m's avatar
frey_m committed
337 338
    PData.setTotalNum(TotalNum);    // set the total atom count
    PData.setLocalNum(LocalNum);    // set the number of local atoms
kraus's avatar
kraus committed
339

frey_m's avatar
AMR:  
frey_m committed
340 341 342 343
    // final check
    if ( LocalNum != LocalNumPerLevel.getLocalNumAllLevel() )
        throw OpalException("BoxLibLayout::update()",
                            "Local #particles disagrees with sum over levels");
kraus's avatar
kraus committed
344

frey_m's avatar
frey_m committed
345
    if ( !PData.isForbidTransform() ) {
346
        // undo domain transformation + undo Lorentz transform
frey_m's avatar
frey_m committed
347
        PData.domainMapping(true);
348
    }
frey_m's avatar
frey_m committed
349 350 351
}


frey_m's avatar
frey_m committed
352
// Function from AMReX adjusted to work with Ippl AmrParticleBase class
frey_m's avatar
frey_m committed
353 354
//get the cell where particle is located - uses AmrParticleBase object and particle id
template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
355 356 357 358
typename BoxLibLayout<T, Dim>::AmrIntVect_t
    BoxLibLayout<T, Dim>::Index(AmrParticleBase< BoxLibLayout<T,Dim> >& p,
                                const unsigned int ip,
                                int lev) const
frey_m's avatar
frey_m committed
359
{
frey_m's avatar
frey_m committed
360
    return Index(p.R[ip], lev);
frey_m's avatar
frey_m committed
361 362 363 364
}

//get the cell where particle is located - uses the particle position vector R
template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
365 366 367
typename BoxLibLayout<T, Dim>::AmrIntVect_t
    BoxLibLayout<T, Dim>::Index(SingleParticlePos_t &R,
                                int lev) const
frey_m's avatar
frey_m committed
368
{
frey_m's avatar
frey_m committed
369 370
    AmrIntVect_t iv;
    const AmrGeometry_t& geom = Geom(lev);
frey_m's avatar
frey_m committed
371 372 373 374 375 376 377 378 379 380

    D_TERM(iv[0]=floor((R[0]-geom.ProbLo(0))/geom.CellSize(0));,
           iv[1]=floor((R[1]-geom.ProbLo(1))/geom.CellSize(1));,
           iv[2]=floor((R[2]-geom.ProbLo(2))/geom.CellSize(2)););

    iv += geom.Domain().smallEnd();

    return iv;
}

381

frey_m's avatar
frey_m committed
382 383 384 385 386 387
template <class T, unsigned Dim>
void BoxLibLayout<T, Dim>::buildLevelMask(int lev, const int ncells) {
    int covered = 0;
    int notcovered = 1;
    int physbnd = 1;
    int interior = 0;
kraus's avatar
kraus committed
388

frey_m's avatar
frey_m committed
389 390
    if ( lev >= (int)masks_m.size() )
        masks_m.resize(lev + 1);
kraus's avatar
kraus committed
391

frey_m's avatar
frey_m committed
392 393
    masks_m[lev].reset(new mask_t(ParticleBoxArray(lev),
                                  ParticleDistributionMap(lev), 1, 1));
kraus's avatar
kraus committed
394

frey_m's avatar
frey_m committed
395
    masks_m[lev]->setVal(1, 1);
kraus's avatar
kraus committed
396

frey_m's avatar
frey_m committed
397 398 399
    mask_t tmp_mask(ParticleBoxArray(lev),
                    ParticleDistributionMap(lev),
                    1, ncells);
kraus's avatar
kraus committed
400

frey_m's avatar
frey_m committed
401
    tmp_mask.setVal(0, ncells);
kraus's avatar
kraus committed
402

frey_m's avatar
frey_m committed
403 404
    tmp_mask.BuildMask(Geom(lev).Domain(), Geom(lev).periodicity(),
                       covered, notcovered, physbnd, interior);
kraus's avatar
kraus committed
405

frey_m's avatar
frey_m committed
406
    tmp_mask.FillBoundary(Geom(lev).periodicity());
kraus's avatar
kraus committed
407

frey_m's avatar
frey_m committed
408 409 410 411
    for (amrex::MFIter mfi(tmp_mask); mfi.isValid(); ++mfi) {
        const AmrBox_t& bx = mfi.validbox();
        const int* lo = bx.loVect();
        const int* hi = bx.hiVect();
kraus's avatar
kraus committed
412

frey_m's avatar
frey_m committed
413 414
        basefab_t& mfab = (*masks_m[lev])[mfi];
        const basefab_t& fab = tmp_mask[mfi];
kraus's avatar
kraus committed
415

frey_m's avatar
frey_m committed
416 417 418 419
        for (int i = lo[0]; i <= hi[0]; ++i) {
            for (int j = lo[1]; j <= hi[1]; ++j) {
                for (int k = lo[2]; k <= hi[2]; ++k) {
                    int total = 0;
kraus's avatar
kraus committed
420

frey_m's avatar
frey_m committed
421 422 423 424 425 426 427 428
                    for (int ii = i - ncells; ii <= i + ncells; ++ii) {
                        for (int jj = j - ncells; jj <= j + ncells; ++jj) {
                            for (int kk = k - ncells; kk <= k + ncells; ++kk) {
                                AmrIntVect_t iv(ii, jj, kk);
                                total += fab(iv);
                            }
                        }
                    }
kraus's avatar
kraus committed
429

frey_m's avatar
frey_m committed
430 431 432 433 434 435 436 437
                    AmrIntVect_t iv(i, j, k);
                    if (total == 0) {
                        mfab(iv) = 0;
                    }
                }
            }
        }
    }
kraus's avatar
kraus committed
438

frey_m's avatar
frey_m committed
439 440 441 442 443
    masks_m[lev]->FillBoundary(Geom(lev).periodicity());
}


template <class T, unsigned Dim>
444
void BoxLibLayout<T, Dim>::clearLevelMask(int lev) {
445
    assert(lev < (int)masks_m.size());
446 447 448 449 450 451
    masks_m[lev].reset(nullptr);
}


template <class T, unsigned Dim>
const std::unique_ptr<typename BoxLibLayout<T, Dim>::mask_t>& 
frey_m's avatar
frey_m committed
452 453 454 455 456 457 458 459
BoxLibLayout<T, Dim>::getLevelMask(int lev) const {
    if ( lev >= (int)masks_m.size() ) {
        throw OpalException("BoxLibLayout::getLevelMask()",
                            "Unable to access level " + std::to_string(lev) + ".");
    }
    return masks_m[lev];
}

frey_m's avatar
frey_m committed
460

frey_m's avatar
frey_m committed
461
// template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
462
// int BoxLibLayout<T, Dim>::getTileIndex(const AmrIntVect_t& iv, const Box& box, Box& tbx) {
frey_m's avatar
frey_m committed
463 464 465 466 467 468 469 470 471 472 473 474 475 476
//     if (do_tiling == false) {
//         tbx = box;
//         return 0;
//     } else {
//         //
//         // This function must be consistent with FabArrayBase::buildTileArray function!!!
//         //
//         auto tiling_1d = [](int i, int lo, int hi, int tilesize,
//                             int& ntile, int& tileidx, int& tlo, int& thi) {
//             int ncells = hi-lo+1;
//             ntile = std::max(ncells/tilesize, 1);
//             int ts_right = ncells/ntile;
//             int ts_left  = ts_right+1;
//             int nleft = ncells - ntile*ts_right;
kraus's avatar
kraus committed
477
//             int ii = i - lo;
frey_m's avatar
frey_m committed
478 479 480 481 482 483 484 485 486 487 488
//             int nbndry = nleft*ts_left;
//             if (ii < nbndry) {
//                 tileidx = ii / ts_left; // tiles on the left of nbndry have size of ts_left
//                 tlo = lo + tileidx * ts_left;
//                 thi = tlo + ts_left - 1;
//             } else {
//                 tileidx = nleft + (ii-nbndry) / ts_right;  // tiles on the right: ts_right
//                 tlo = lo + tileidx * ts_right + nleft;
//                 thi = tlo + ts_right - 1;
//             }
//         };
frey_m's avatar
frey_m committed
489 490 491
//         const AmrIntVect_t& small = box.smallEnd();
//         const AmrIntVect_t& big   = box.bigEnd();
//         AmrIntVect_t ntiles, ivIndex, tilelo, tilehi;
kraus's avatar
kraus committed
492
//
frey_m's avatar
frey_m committed
493 494 495
//         D_TERM(int iv0 = std::min(std::max(iv[0], small[0]), big[0]);,
//                int iv1 = std::min(std::max(iv[1], small[1]), big[1]);,
//                int iv2 = std::min(std::max(iv[2], small[2]), big[2]););
kraus's avatar
kraus committed
496
//
frey_m's avatar
frey_m committed
497 498 499
//         D_TERM(tiling_1d(iv0, small[0], big[0], tile_size[0], ntiles[0], ivIndex[0], tilelo[0], tilehi[0]);,
//                tiling_1d(iv1, small[1], big[1], tile_size[1], ntiles[1], ivIndex[1], tilelo[1], tilehi[1]);,
//                tiling_1d(iv2, small[2], big[2], tile_size[2], ntiles[2], ivIndex[2], tilelo[2], tilehi[2]););
kraus's avatar
kraus committed
500
//
frey_m's avatar
frey_m committed
501
//         tbx = Box(tilelo, tilehi);
kraus's avatar
kraus committed
502
//
frey_m's avatar
frey_m committed
503 504 505 506 507 508
//         return D_TERM(ivIndex[0], + ntiles[0]*ivIndex[1], + ntiles[0]*ntiles[1]*ivIndex[2]);
//     }
// }


//sets the grid and level where particle belongs - returns false if particle is outside the domain
frey_m's avatar
frey_m committed
509
template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
510 511 512 513 514
bool BoxLibLayout<T, Dim>::Where(AmrParticleBase< BoxLibLayout<T,Dim> >& p,
                                 const unsigned int ip,
                                 int lev_min,
                                 int lev_max,
                                 int nGrow) const
frey_m's avatar
frey_m committed
515 516
{

frey_m's avatar
frey_m committed
517 518
    if (lev_max == -1)
        lev_max = finestLevel();
kraus's avatar
kraus committed
519

frey_m's avatar
frey_m committed
520
    BL_ASSERT(lev_max <= finestLevel());
kraus's avatar
kraus committed
521

frey_m's avatar
frey_m committed
522
    BL_ASSERT(nGrow == 0 || (nGrow >= 0 && lev_min == lev_max));
frey_m's avatar
frey_m committed
523

frey_m's avatar
frey_m committed
524
    std::vector< std::pair<int, AmrBox_t> > isects;
frey_m's avatar
frey_m committed
525

frey_m's avatar
frey_m committed
526
    for (int lev = lev_max; lev >= lev_min; lev--)
frey_m's avatar
frey_m committed
527
    {
frey_m's avatar
frey_m committed
528 529
        const AmrIntVect_t& iv = Index(p, ip, lev);
        const AmrGrid_t& ba = ParticleBoxArray(lev);
frey_m's avatar
frey_m committed
530 531
        BL_ASSERT(ba.ixType().cellCentered());

kraus's avatar
kraus committed
532
        if (lev == (int)p.Level[ip]) {
frey_m's avatar
frey_m committed
533
            // The fact that we are here means this particle does not belong to any finer grids.
frey_m's avatar
frey_m committed
534 535 536
            if (0 <= p.Grid[ip] && p.Grid[ip] < ba.size())
            {
                const AmrBox_t& bx = ba.getCellCenteredBox(p.Grid[ip]);
frey_m's avatar
frey_m committed
537
                const AmrBox_t& gbx = amrex::grow(bx,nGrow);
frey_m's avatar
frey_m committed
538 539 540 541 542 543 544 545
                if (gbx.contains(iv))
                {
//                     if (bx != pld.m_gridbox || !pld.m_tilebox.contains(iv)) {
//                         pld.m_tile = getTileIndex(iv, bx, pld.m_tilebox);
//                         pld.m_gridbox = bx;
//                     }
                    return true;
                }
frey_m's avatar
frey_m committed
546 547
            }
        }
kraus's avatar
kraus committed
548

frey_m's avatar
frey_m committed
549
        ba.intersections(AmrBox_t(iv, iv), isects, true, nGrow);
kraus's avatar
kraus committed
550

frey_m's avatar
frey_m committed
551 552
        if (!isects.empty())
        {
553 554
            p.Level[ip]  = lev;
            p.Grid[ip] = isects[0].first;
frey_m's avatar
frey_m committed
555 556 557 558 559 560 561 562

            return true;
        }
    }
    return false;
}


frey_m's avatar
frey_m committed
563
//Function from AMReX adjusted to work with Ippl AmrParticleBase class
frey_m's avatar
frey_m committed
564 565 566
//Checks/sets whether the particle has crossed a periodic boundary in such a way
//that it is on levels lev_min and higher.
template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
567 568 569 570
bool BoxLibLayout<T, Dim>::EnforcePeriodicWhere (AmrParticleBase< BoxLibLayout<T,Dim> >& p,
                                                      const unsigned int ip,
                                                      int lev_min,
                                                      int lev_max) const
frey_m's avatar
frey_m committed
571
{
frey_m's avatar
frey_m committed
572
    if (!Geom(0).isAnyPeriodic()) return false;
frey_m's avatar
frey_m committed
573

frey_m's avatar
frey_m committed
574 575
    if (lev_max == -1)
        lev_max = finestLevel();
frey_m's avatar
frey_m committed
576

frey_m's avatar
frey_m committed
577
    BL_ASSERT(lev_max <= finestLevel());
frey_m's avatar
frey_m committed
578 579 580 581 582
    //
    // Create a copy "dummy" particle to check for periodic outs.
    //
    SingleParticlePos_t R = p.R[ip];

frey_m's avatar
frey_m committed
583
    if (PeriodicShift(R))
frey_m's avatar
frey_m committed
584
    {
frey_m's avatar
frey_m committed
585
        std::vector< std::pair<int, AmrBox_t> > isects;
frey_m's avatar
frey_m committed
586

frey_m's avatar
frey_m committed
587
        for (int lev = lev_max; lev >= lev_min; lev--)
frey_m's avatar
frey_m committed
588
        {
frey_m's avatar
frey_m committed
589 590
            const AmrIntVect_t& iv = Index(R, lev);
            const AmrGrid_t& ba = ParticleBoxArray(lev);
kraus's avatar
kraus committed
591

frey_m's avatar
frey_m committed
592
            ba.intersections(AmrBox_t(iv,iv),isects,true,0);
frey_m's avatar
frey_m committed
593 594 595 596 597 598 599

            if (!isects.empty())
            {
                D_TERM(p.R[ip][0] = R[0];,
                       p.R[ip][1] = R[1];,
                       p.R[ip][2] = R[2];);

600 601
                p.Level[ip]  = lev;
                p.Grid[ip] = isects[0].first;
frey_m's avatar
frey_m committed
602 603 604 605 606 607 608 609 610 611

                return true;
            }
        }
    }

    return false;
}


frey_m's avatar
frey_m committed
612
// Function from AMReX adjusted to work with Ippl AmrParticleBase class
frey_m's avatar
frey_m committed
613 614
// Returns true if the particle was shifted.
template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
615
bool BoxLibLayout<T, Dim>::PeriodicShift (SingleParticlePos_t R) const
frey_m's avatar
frey_m committed
616
{
frey_m's avatar
frey_m committed
617
  //
frey_m's avatar
frey_m committed
618 619 620 621 622
    // This routine should only be called when Where() returns false.
    //
    //
    // We'll use level 0 stuff since ProbLo/ProbHi are the same for every level.
    //
frey_m's avatar
frey_m committed
623 624 625
    const AmrGeometry_t& geom    = Geom(0);
    const AmrBox_t&      dmn     = geom.Domain();
    const AmrIntVect_t&  iv      = Index(R, 0);
kraus's avatar
kraus committed
626
    bool            shifted = false;
frey_m's avatar
frey_m committed
627

628
    for (int i = 0; i < AMREX_SPACEDIM; i++)
frey_m's avatar
frey_m committed
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683
    {
        if (!geom.isPeriodic(i)) continue;

        if (iv[i] > dmn.bigEnd(i))
        {
            if (R[i] == geom.ProbHi(i))
                //
                // Don't let particles lie exactly on the domain face.
                // Force the particle to be outside the domain so the
                // periodic shift will bring it back inside.
                //
                R[i] += .125*geom.CellSize(i);

            R[i] -= geom.ProbLength(i);

            if (R[i] <= geom.ProbLo(i))
                //
                // This can happen due to precision issues.
                //
                R[i] += .125*geom.CellSize(i);

            BL_ASSERT(R[i] >= geom.ProbLo(i));

            shifted = true;
        }
        else if (iv[i] < dmn.smallEnd(i))
        {
            if (R[i] == geom.ProbLo(i))
                //
                // Don't let particles lie exactly on the domain face.
                // Force the particle to be outside the domain so the
                // periodic shift will bring it back inside.
                //
                R[i] -= .125*geom.CellSize(i);

            R[i] += geom.ProbLength(i);

            if (R[i] >= geom.ProbHi(i))
                //
                // This can happen due to precision issues.
                //
                R[i] -= .125*geom.CellSize(i);

            BL_ASSERT(R[i] <= geom.ProbHi(i));

            shifted = true;
        }
    }
    //
    // The particle may still be outside the domain in the case
    // where we aren't periodic on the face out which it travelled.
    //
    return shifted;
}

frey_m's avatar
frey_m committed
684

685 686 687 688
template <class T, unsigned Dim>
void BoxLibLayout<T, Dim>::locateParticle(
    AmrParticleBase< BoxLibLayout<T,Dim> >& p,
    const unsigned int ip,
frey_m's avatar
frey_m committed
689
    int lev_min, int lev_max, int nGrow) const
690 691 692 693 694 695 696
{
    bool outside = D_TERM( p.R[ip](0) <  AmrGeometry_t::ProbLo(0)
                        || p.R[ip](0) >= AmrGeometry_t::ProbHi(0),
                        || p.R[ip](1) <  AmrGeometry_t::ProbLo(1)
                        || p.R[ip](1) >= AmrGeometry_t::ProbHi(1),
                        || p.R[ip](2) <  AmrGeometry_t::ProbLo(2)
                        || p.R[ip](2) >= AmrGeometry_t::ProbHi(2));
kraus's avatar
kraus committed
697

frey_m's avatar
frey_m committed
698
    bool success = false;
kraus's avatar
kraus committed
699

700 701 702 703 704 705 706 707 708
    if (outside)
    {
        // Note that EnforcePeriodicWhere may shift the particle if it is successful.
        success = EnforcePeriodicWhere(p, ip, lev_min, lev_max);
        if (!success && lev_min == 0)
        {
            // The particle has left the domain; invalidate it.
            p.destroy(1, ip);
            success = true;
kraus's avatar
kraus committed
709

710 711 712 713 714
            /* We shouldn't lose particles since they are mapped to be within
             * [-1, 1]^3.
             */
            throw OpalException("BoxLibLayout::locateParticle()",
                                "We're losing particles although we shouldn't");
kraus's avatar
kraus committed
715

716 717 718 719 720 721
        }
    }
    else
    {
        success = Where(p, ip, lev_min, lev_max);
    }
kraus's avatar
kraus committed
722

723 724 725 726
    if (!success)
    {
        success = (nGrow > 0) && Where(p, ip, lev_min, lev_min, nGrow);
    }
kraus's avatar
kraus committed
727

728 729
    if (!success)
    {
frey_m's avatar
frey_m committed
730
        std::stringstream ss;
frey_m's avatar
frey_m committed
731
        ss << "Invalid particle with ID " << ip << " at position " << p.R[ip] << ".";
frey_m's avatar
frey_m committed
732
        throw OpalException("BoxLibLayout::locateParticle()", ss.str());
733 734
    }
}
frey_m's avatar
frey_m committed
735

frey_m's avatar
frey_m committed
736

frey_m's avatar
frey_m committed
737 738 739 740 741 742 743 744 745 746 747 748 749
// overwritten functions
template <class T, unsigned Dim>
bool BoxLibLayout<T, Dim>::LevelDefined (int level) const {
    return level <= this->maxLevel_m && !m_ba[level].empty() && !m_dmap[level].empty();
}


template <class T, unsigned Dim>
int BoxLibLayout<T, Dim>::finestLevel () const {
    return this->finestLevel_m;
}


frey_m's avatar
frey_m committed
750
template <class T, unsigned Dim>
frey_m's avatar
frey_m committed
751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766
int BoxLibLayout<T, Dim>::maxLevel () const {
    return this->maxLevel_m;
}


template <class T, unsigned Dim>
typename BoxLibLayout<T, Dim>::AmrIntVect_t
    BoxLibLayout<T, Dim>::refRatio (int level) const
{
    return refRatio_m[level];
}


template <class T, unsigned Dim>
int BoxLibLayout<T, Dim>::MaxRefRatio (int level) const {
    int maxval = 0;
kraus's avatar
kraus committed
767
    for (int n = 0; n<AMREX_SPACEDIM; n++)
frey_m's avatar
frey_m committed
768 769 770 771 772 773 774 775 776
        maxval = std::max(maxval, refRatio_m[level][n]);
    return maxval;
}


template <class T, unsigned Dim>
void BoxLibLayout<T, Dim>::initBaseBox_m(int nGridPoints,
                                         int maxGridSize,
                                         double dh)
frey_m's avatar
frey_m committed
777
{
frey_m's avatar
frey_m committed
778
    // physical box (in meters)
frey_m's avatar
frey_m committed
779
    AmrDomain_t real_box;
780
    for (int d = 0; d < AMREX_SPACEDIM; ++d) {
781 782 783 784 785 786

        assert(lowerBound[d] < 0);
        assert(upperBound[d] > 0);

        real_box.setLo(d, lowerBound[d] * (1.0 + dh));
        real_box.setHi(d, upperBound[d] * (1.0 + dh));
frey_m's avatar
frey_m committed
787
    }
kraus's avatar
kraus committed
788

frey_m's avatar
frey_m committed
789
    AmrGeometry_t::ProbDomain(real_box);
kraus's avatar
kraus committed
790

frey_m's avatar
frey_m committed
791
    // define underlying box for physical domain
kraus's avatar
kraus committed
792 793
    AmrIntVect_t domain_lo(0 , 0, 0);
    AmrIntVect_t domain_hi(nGridPoints - 1, nGridPoints - 1, nGridPoints - 1);
frey_m's avatar
frey_m committed
794
    const AmrBox_t domain(domain_lo, domain_hi);
frey_m's avatar
frey_m committed
795 796 797 798

    // use Cartesian coordinates
    int coord = 0;

frey_m's avatar
frey_m committed
799
    // Dirichlet boundary conditions
800
    int is_per[AMREX_SPACEDIM];
kraus's avatar
kraus committed
801
    for (int i = 0; i < AMREX_SPACEDIM; i++)
frey_m's avatar
frey_m committed
802
        is_per[i] = 0;
kraus's avatar
kraus committed
803

frey_m's avatar
frey_m committed
804
    AmrGeometry_t geom;
frey_m's avatar
frey_m committed
805 806
    geom.define(domain, &real_box, coord, is_per);

frey_m's avatar
frey_m committed
807
    AmrGrid_t ba;
frey_m's avatar
frey_m committed
808 809 810
    ba.define(domain);
    // break the BoxArrays at both levels into max_grid_size^3 boxes
    ba.maxSize(maxGridSize);
kraus's avatar
kraus committed
811

frey_m's avatar
frey_m committed
812
    AmrProcMap_t dmap;
frey_m's avatar
frey_m committed
813
    dmap.define(ba, Ippl::getNodes());
kraus's avatar
kraus committed
814

frey_m's avatar
frey_m committed
815 816 817
    // set protected ParGDB member variables
    this->m_geom.resize(1);
    this->m_geom[0] = geom;
kraus's avatar
kraus committed
818

frey_m's avatar
frey_m committed
819 820
    this->m_dmap.resize(1);
    this->m_dmap[0] = dmap;
kraus's avatar
kraus committed
821

frey_m's avatar
frey_m committed
822 823
    this->m_ba.resize(1);
    this->m_ba[0] = ba;
kraus's avatar
kraus committed
824

frey_m's avatar
frey_m committed
825 826 827
    this->m_nlevels = ba.size();
}

kraus's avatar
kraus committed
828
#endif