AmrBoxLib.cpp 50.5 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
//
// Class AmrBoxLib
//   Concrete AMR object. It is based on the AMReX library
//   (cf. https://amrex-codes.github.io/ or https://ccse.lbl.gov/AMReX/).
//   AMReX is the successor of BoxLib. This class represents the interface
//   to AMReX and the AMR framework in OPAL. It implements the functions of
//   the AmrObject class.
//
// Copyright (c) 2016 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
25 26
#include "AmrBoxLib.h"

27
#include "Algorithms/AmrPartBunch.h"
frey_m's avatar
AMR:  
frey_m committed
28
#include "Structure/FieldSolver.h"
frey_m's avatar
frey_m committed
29
#include "Solvers/PoissonSolver.h"
snuverink_j's avatar
snuverink_j committed
30
#include "Utility/PAssert.h"
31

32
#include "Amr/AmrYtWriter.h"
snuverink_j's avatar
snuverink_j committed
33

frey_m's avatar
frey_m committed
34
#include <AMReX_MultiFabUtil.H>
35

frey_m's avatar
frey_m committed
36
#include <AMReX_ParmParse.H> // used in initialize function
frey_m's avatar
frey_m committed
37
#include <AMReX_BCUtil.H>
38

frey_m's avatar
frey_m committed
39 40
#include <map>

41 42
extern Inform* gmsg;

frey_m's avatar
frey_m committed
43 44 45

AmrBoxLib::AmrBoxLib(const AmrDomain_t& domain,
                     const AmrIntArray_t& nGridPts,
frey_m's avatar
frey_m committed
46
                     int maxLevel,
frey_m's avatar
frey_m committed
47
                     AmrPartBunch* bunch_p)
frey_m's avatar
frey_m committed
48 49 50 51 52 53 54 55 56 57
    : AmrObject()
    , amrex::AmrMesh(&domain, maxLevel, nGridPts, 0 /* cartesian */)
    , bunch_mp(bunch_p)
    , layout_mp(static_cast<AmrLayout_t*>(&bunch_p->getLayout()))
    , rho_m(maxLevel + 1)
    , phi_m(maxLevel + 1)
    , efield_m(maxLevel + 1)
    , meshScaling_m(Vector_t(1.0, 1.0, 1.0))
    , isFirstTagging_m(maxLevel + 1, true)
    , isPoissonSolved_m(false)
58
{
frey_m's avatar
frey_m committed
59 60 61 62 63 64 65 66 67
    /*
     * The layout needs to know how many levels we can make.
     */
    layout_mp->resize(maxLevel);
    
    initBaseLevel_m(nGridPts);
    
    // set mesh spacing of bunch
    updateMesh();
68 69 70
}


frey_m's avatar
frey_m committed
71
std::unique_ptr<AmrBoxLib> AmrBoxLib::create(const AmrInfo& info,
frey_m's avatar
AMR:  
frey_m committed
72 73 74 75 76 77 78 79 80 81 82
                                             AmrPartBunch* bunch_p)
{
    /* The bunch is initialized first with a Geometry,
     * BoxArray and DistributionMapping on
     * the base level (level = 0). Thus, we take the domain specified there in
     * order to create the Amr object.
     */
    AmrLayout_t* layout_p = static_cast<AmrLayout_t*>(&bunch_p->getLayout());
    AmrDomain_t domain = layout_p->Geom(0).ProbDomain();
    
    AmrIntArray_t nGridPts = {
frey_m's avatar
frey_m committed
83 84 85
        info.grid[0],
        info.grid[1],
        info.grid[2]
frey_m's avatar
AMR:  
frey_m committed
86 87 88 89 90 91 92
    };
    
    int maxlevel = info.maxlevel;
    
    /*
     * further attributes are given by the BoxLib's ParmParse class.
     */
frey_m's avatar
frey_m committed
93
    initParmParse_m(info, layout_p);
frey_m's avatar
AMR:  
frey_m committed
94 95 96 97 98 99 100 101 102
    
    return std::unique_ptr<AmrBoxLib>(new AmrBoxLib(domain,
                                                    nGridPts,
                                                    maxlevel,
                                                    bunch_p
                                                   )
                                     );
}

103

104 105 106 107 108 109 110
void AmrBoxLib::regrid(double time) {
    IpplTimings::startTimer(this->amrRegridTimer_m);

    *gmsg << "* Start regriding:" << endl
          << "*     Old finest level: "
          << finest_level << endl;

frey_m's avatar
frey_m committed
111 112
    this->preRegrid_m();

113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
    /* ATTENTION: The bunch might be updated during
     * the regrid process!
     * We regrid from base level 0 up to the finest level.
     */
    int old_finest = finest_level;
    
    int lev_top = std::min(finest_level, max_level - 1);
    for (int i = 0; i <= lev_top; ++i) {
        this->doRegrid_m(i, time);
        lev_top = std::min(finest_level, max_level - 1);
    }

    this->postRegrid_m(old_finest);

    *gmsg << "*     New finest level: "
          << finest_level << endl
          << "* Finished regriding" << endl;

    IpplTimings::stopTimer(this->amrRegridTimer_m);
}


void AmrBoxLib::getGridStatistics(std::map<int, long>& gridPtsPerCore,
136 137 138 139
                                  std::vector<int>& gridsPerLevel) const
{
    typedef std::vector<int> container_t;
    
140
    gridPtsPerCore.clear();
141 142 143 144 145 146 147 148
    gridsPerLevel.clear();
    
    gridsPerLevel.resize(max_level + 1);
    
    for (int lev = 0; lev <= finest_level; ++lev) {
        /* container index: box
         * container value: cores that owns box
         */
frey_m's avatar
frey_m committed
149
        const container_t& pmap = this->dmap[lev].ProcessorMap();
150
        const AmrGrid_t& ba = this->grids[lev];
151 152
        
        gridsPerLevel[lev] = pmap.size();
153 154 155 156 157

        // iterate over all boxes
        for (unsigned int i = 0; i < ba.size(); ++i) {
            gridPtsPerCore[pmap[i]] += ba[i].numPts();
        }
158 159 160 161
    }
}


frey_m's avatar
frey_m committed
162
void AmrBoxLib::initFineLevels() {
frey_m's avatar
frey_m committed
163
    if ( bunch_mp->getNumBunch() > 1 && !refined_m ) {
frey_m's avatar
frey_m committed
164 165 166 167 168 169 170 171 172 173 174 175 176 177
        *gmsg << "* Initialization of all levels" << endl;
        
        AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase();
        
        /* we do an explicit domain mapping of the particles and then
         * forbid it during the regrid process, this way it's only
         * executed ones --> saves computation
         */
        bool isForbidTransform = amrpbase_p->isForbidTransform();
        
        if ( !isForbidTransform ) {
            amrpbase_p->domainMapping();
            amrpbase_p->setForbidTransform(true);
        }
178

frey_m's avatar
frey_m committed
179
        if ( max_level > 0 ) {
180
            this->regrid(bunch_mp->getT() * 1.0e9 /*time [ns] */);
frey_m's avatar
frey_m committed
181
        }
frey_m's avatar
frey_m committed
182

frey_m's avatar
frey_m committed
183 184 185 186 187 188 189 190 191 192 193
        if ( !isForbidTransform ) {
            amrpbase_p->setForbidTransform(false);
            // map particles back
            amrpbase_p->domainMapping(true);
        }
        
        *gmsg << "* Initialization done." << endl;
        
        refined_m = true;
    }
}
frey_m's avatar
AMR:  
frey_m committed
194

195

196 197 198 199 200
AmrBoxLib::VectorPair_t AmrBoxLib::getEExtrema() {
    Vector_t maxE = Vector_t(0, 0, 0), minE = Vector_t(0, 0, 0);
    /* check for maximum / minimum values over all levels
     * and all components, i.e. Ex, Ey, Ez
     */
frey_m's avatar
frey_m committed
201
    for (unsigned int lev = 0; lev < efield_m.size(); ++lev) {
202 203 204 205 206
        for (std::size_t i = 0; i < bunch_mp->Dimension; ++i) {
            /* calls:
             *  - max(comp, nghost (to search), local)
             *  - min(cmop, nghost, local)
             */
207
            double max = efield_m[lev][i]->max(0, false);
208 209
            maxE[i] = (maxE[i] < max) ? max : maxE[i];
            
210
            double min = efield_m[lev][i]->min(0, false);
211 212 213 214 215 216 217 218
            minE[i] = (minE[i] > min) ? min : minE[i];
        }
    }
    
    return VectorPair_t(maxE, minE);
}


frey_m's avatar
frey_m committed
219
double AmrBoxLib::getRho(int /*x*/, int /*y*/, int /*z*/) {
220
    //TODO
221
    throw OpalException("AmrBoxLib::getRho(x, y, z)", "Not yet Implemented.");
222 223 224 225 226
    return 0.0;
}


void AmrBoxLib::computeSelfFields() {
227 228
    //TODO
    throw OpalException("AmrBoxLib::computeSelfFields", "Not yet Implemented.");
229 230 231
}


frey_m's avatar
frey_m committed
232
void AmrBoxLib::computeSelfFields(int /*bin*/) {
233
    //TODO
234
    throw OpalException("AmrBoxLib::computeSelfFields(int bin)", "Not yet Implemented.");
235 236 237
}


frey_m's avatar
frey_m committed
238 239 240 241 242
void AmrBoxLib::computeSelfFields_cycl(double gamma) {   
    /*
     * The potential is not scaled according to domain modification.
     * 
     */
frey_m's avatar
frey_m committed
243 244
    if ( !bunch_mp->hasFieldSolver() )
        return;
frey_m's avatar
frey_m committed
245 246 247 248
    
    /*
     * scatter charges onto grid
     */
249
    AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase();
250
    
frey_m's avatar
frey_m committed
251 252
    /// Lorentz transformation
    /// In particle rest frame, the longitudinal length (y for cyclotron) enlarged
253
    bunch_mp->updateLorentzFactor(gamma);
254
    
255
    // map on Amr domain + Lorentz transform
frey_m's avatar
frey_m committed
256
    double scalefactor = amrpbase_p->domainMapping();
frey_m's avatar
frey_m committed
257
    
258 259 260 261
    amrpbase_p->setForbidTransform(true);
    amrpbase_p->update();
    amrpbase_p->setForbidTransform(false);
    
frey_m's avatar
frey_m committed
262
    double invGamma = 1.0 / gamma;
frey_m's avatar
frey_m committed
263
    int nLevel = finest_level + 1;
frey_m's avatar
frey_m committed
264

frey_m's avatar
frey_m committed
265
    double l0norm = this->solvePoisson_m();
frey_m's avatar
frey_m committed
266

267 268 269 270 271
    /* apply scale of electric-field in order to undo the transformation
     * + undo normalization
     */
    for (int i = 0; i <= finestLevel(); ++i) {
        this->phi_m[i]->mult(scalefactor * l0norm, 0, 1);
272 273 274
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            this->efield_m[i][j]->mult(scalefactor * scalefactor * l0norm, 0, 1);
        }
275
    }
frey_m's avatar
frey_m committed
276

frey_m's avatar
frey_m committed
277
    for (int i = 0; i <= finest_level; ++i) {
278 279 280 281 282
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            if ( this->efield_m[i][j]->contains_nan(false) )
                throw OpalException("AmrBoxLib::computeSelfFields_cycl(double gamma) ",
                                    "Ef: NANs at level " + std::to_string(i) + ".");
        }
frey_m's avatar
frey_m committed
283
    }
frey_m's avatar
frey_m committed
284

frey_m's avatar
frey_m committed
285
    amrpbase_p->gather(bunch_mp->Ef, this->efield_m, bunch_mp->R, 0, finest_level);
286
    
287
    // undo domain change + undo Lorentz transform
frey_m's avatar
frey_m committed
288
    amrpbase_p->domainMapping(true);
289 290 291
    
    /// Back Lorentz transformation
    bunch_mp->Ef *= Vector_t(gamma,
frey_m's avatar
frey_m committed
292
                             1.0,
293 294 295 296 297 298 299 300 301 302 303 304 305 306 307
                             gamma);
    
    /// calculate coefficient
    // Relativistic E&M says gamma*v/c^2 = gamma*beta/c = sqrt(gamma*gamma-1)/c
    // but because we already transformed E_trans into the moving frame we have to
    // add 1/gamma so we are using the E_trans from the rest frame -DW
    double betaC = std::sqrt(gamma * gamma - 1.0) * invGamma / Physics::c;
    
    /// calculate B field from E field
    bunch_mp->Bf(0) =  betaC * bunch_mp->Ef(2);
    bunch_mp->Bf(2) = -betaC * bunch_mp->Ef(0);
    
    /*
     * dumping only
     */
frey_m's avatar
frey_m committed
308 309 310
    
    if ( !(bunch_mp->getLocalTrackStep()  % Options::amrYtDumpFreq) ) {
        AmrYtWriter ytWriter(bunch_mp->getLocalTrackStep());
311
    
312 313
        AmrIntArray_t rr(finest_level);
        for (int i = 0; i < finest_level; ++i)
314 315 316 317 318 319 320 321 322
            rr[i] = this->MaxRefRatio(i);
        
        double time = bunch_mp->getT(); // in seconds
        
        // we need to undo coefficient when writing charge density
        for (int i = 0; i <= finest_level; ++i)
            this->rho_m[i]->mult(- Physics::epsilon_0 * l0norm, 0, 1);
        
        
frey_m's avatar
frey_m committed
323 324
        ytWriter.writeFields(rho_m, phi_m, efield_m, rr, this->geom,
                             nLevel, time, scalefactor);
325
        
326
        ytWriter.writeBunch(bunch_mp, time, gamma);
327
    }
328 329 330
}


331
void AmrBoxLib::computeSelfFields_cycl(int bin) {
frey_m's avatar
frey_m committed
332

frey_m's avatar
AMR:  
frey_m committed
333 334
    if ( !bunch_mp->hasFieldSolver() )
        return;
frey_m's avatar
frey_m committed
335
    
frey_m's avatar
AMR:  
frey_m committed
336 337 338
    /// get gamma of this bin
    double gamma = bunch_mp->getBinGamma(bin);
    
339 340 341 342 343 344
    /* relativistic factor is always gamma >= 1
     * --> if no particle --> gamma = 0 --> leave computation
     */
    if ( gamma < 1.0 )
        return;

frey_m's avatar
AMR:  
frey_m committed
345 346 347 348
    /*
     * scatter charges onto grid
     */
    AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase();
349
    
frey_m's avatar
frey_m committed
350
    // Lorentz transformation: apply Lorentz factor of that particle bin
351
    bunch_mp->updateLorentzFactor(bin);
352
    
353
    // map on Amr domain + Lorentz transform 
frey_m's avatar
AMR:  
frey_m committed
354 355
    double scalefactor = amrpbase_p->domainMapping();
    
356 357 358
    amrpbase_p->setForbidTransform(true);
    
    amrpbase_p->update();
frey_m's avatar
frey_m committed
359 360

    if ( !(bunch_mp->getLocalTrackStep() % Options::amrRegridFreq) ) {
361
        this->regrid(bunch_mp->getT() * 1.0e9 /*time [ns] */);
frey_m's avatar
frey_m committed
362 363
    }

364 365
    amrpbase_p->setForbidTransform(false);
    
frey_m's avatar
AMR:  
frey_m committed
366 367
    /// scatter particles charge onto grid.
    /// from charge (C) to charge density (C/m^3).
368
    amrpbase_p->scatter(bunch_mp->Q, this->rho_m, bunch_mp->R, 0, finest_level, bunch_mp->Bin, bin);
frey_m's avatar
AMR:  
frey_m committed
369 370 371 372 373
    
    /// Lorentz transformation
    // In particle rest frame, the longitudinal length (y for cyclotron) enlarged
    // calculate Possion equation (with coefficient: -1/(eps))
    for (int i = 0; i <= finest_level; ++i) {
frey_m's avatar
frey_m committed
374
        this->rho_m[i]->mult(-1.0 / (Physics::epsilon_0), 0 /*comp*/, 1 /*ncomp*/);
frey_m's avatar
frey_m committed
375
        
frey_m's avatar
AMR:  
frey_m committed
376 377 378 379 380
        if ( this->rho_m[i]->contains_nan(false) )
            throw OpalException("AmrBoxLib::computeSelfFields_cycl(int bin) ",
                                "NANs at level " + std::to_string(i) + ".");
    }
    
381 382 383 384 385 386 387 388 389
    // find maximum and normalize each level (faster convergence)
    double l0norm = 1.0;
    for (int i = 0; i <= finest_level; ++i)
        l0norm = std::max(l0norm, this->rho_m[i]->norm0(0));
    
    for (int i = 0; i <= finest_level; ++i) {
        this->rho_m[i]->mult(1.0 / l0norm, 0, 1);
    }
    
frey_m's avatar
frey_m committed
390
    // mesh scaling for solver
frey_m's avatar
frey_m committed
391
    meshScaling_m = Vector_t(1.0, 1.0, 1.0);
frey_m's avatar
frey_m committed
392
    
frey_m's avatar
AMR:  
frey_m committed
393
    PoissonSolver *solver = bunch_mp->getFieldSolver();
frey_m's avatar
frey_m committed
394

frey_m's avatar
AMR:  
frey_m committed
395 396
    IpplTimings::startTimer(this->amrSolveTimer_m);
    
397 398 399 400 401 402 403
    for (int i = 0; i <= finest_level; ++i) {
        phi_m[i]->setVal(0.0, 0, phi_m[i]->nComp(), phi_m[i]->nGrow());
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            efield_m[i][j]->setVal(0.0, 0, efield_m[i][j]->nComp(), efield_m[i][j]->nGrow());
        }
    }
    
frey_m's avatar
frey_m committed
404 405
    // in case of binning we reset phi every time
    solver->solve(rho_m, phi_m, efield_m, 0, finest_level, false);
frey_m's avatar
AMR:  
frey_m committed
406 407 408
    
    IpplTimings::stopTimer(this->amrSolveTimer_m);
    
409 410 411
    // make sure ghost cells are filled
    for (int i = 0; i <= finest_level; ++i) {
        phi_m[i]->FillBoundary(geom[i].periodicity());
412 413 414
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            efield_m[i][j]->FillBoundary(geom[i].periodicity());
        }
415 416
    }
    
frey_m's avatar
frey_m committed
417
    this->fillPhysbc_m(*(this->phi_m[0]), 0);
frey_m's avatar
frey_m committed
418
    
frey_m's avatar
AMR:  
frey_m committed
419
    
420 421 422 423 424
    /* apply scale of electric-field in order to undo the transformation
     * + undo normalization
     */
    for (int i = 0; i <= finestLevel(); ++i) {
        this->phi_m[i]->mult(scalefactor * l0norm, 0, 1);
425 426 427
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            this->efield_m[i][j]->mult(scalefactor * scalefactor * l0norm, 0, 1);
        }
428
    }
frey_m's avatar
AMR:  
frey_m committed
429 430

    for (int i = 0; i <= finest_level; ++i) {
431 432 433 434 435 436 437
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            if ( this->efield_m[i][j]->contains_nan(false) )
                throw OpalException("AmrBoxLib::computeSelfFields_cycl(double gamma) ",
                                    "Ef: NANs at level " + std::to_string(i) + ".");
        }
    }
    
frey_m's avatar
AMR:  
frey_m committed
438 439
    amrpbase_p->gather(bunch_mp->Eftmp, this->efield_m, bunch_mp->R, 0, finest_level);
    
440
    // undo domain change + undo Lorentz transform
frey_m's avatar
AMR:  
frey_m committed
441 442 443
    amrpbase_p->domainMapping(true);
    
    /// Back Lorentz transformation
frey_m's avatar
frey_m committed
444
    bunch_mp->Eftmp *= Vector_t(gamma, 1.0, gamma);
frey_m's avatar
frey_m committed
445

frey_m's avatar
AMR:  
frey_m committed
446
    /// Calculate coefficient
frey_m's avatar
frey_m committed
447
    double betaC = std::sqrt(gamma * gamma - 1.0) / gamma / Physics::c;
frey_m's avatar
frey_m committed
448

frey_m's avatar
AMR:  
frey_m committed
449 450 451
    /// Calculate B_bin field from E_bin field accumulate B and E field
    bunch_mp->Bf(0) += betaC * bunch_mp->Eftmp(2);
    bunch_mp->Bf(2) -= betaC * bunch_mp->Eftmp(0);
frey_m's avatar
frey_m committed
452

frey_m's avatar
AMR:  
frey_m committed
453
    bunch_mp->Ef += bunch_mp->Eftmp;
frey_m's avatar
frey_m committed
454

frey_m's avatar
AMR:  
frey_m committed
455 456 457
    /*
     * dumping only
     */
frey_m's avatar
frey_m committed
458
    if ( !(bunch_mp->getLocalTrackStep()  % Options::amrYtDumpFreq) ) {
frey_m's avatar
AMR:  
frey_m committed
459
        AmrYtWriter ytWriter(bunch_mp->getLocalTrackStep(), bin);
frey_m's avatar
frey_m committed
460 461 462
        
        int nLevel = finest_level + 1;
        
463
        AmrIntArray_t rr(finest_level);
464 465 466 467 468 469 470 471 472 473
        for (int i = 0; i < finest_level; ++i)
            rr[i] = this->MaxRefRatio(i);
        
        double time = bunch_mp->getT(); // in seconds
        
        // we need to undo coefficient when writing charge density
        for (int i = 0; i <= finest_level; ++i)
            this->rho_m[i]->mult(- Physics::epsilon_0 * l0norm, 0, 1);
        
        
frey_m's avatar
frey_m committed
474 475
        ytWriter.writeFields(rho_m, phi_m, efield_m, rr, this->geom,
                             nLevel, time, scalefactor);
476
        
477
        ytWriter.writeBunch(bunch_mp, time, gamma);
478
    }
frey_m's avatar
frey_m committed
479 480

    isPoissonSolved_m = true;
481 482 483 484
}


void AmrBoxLib::updateMesh() {
frey_m's avatar
frey_m committed
485 486
    //FIXME What about resizing mesh, i.e. geometry?    
    const AmrReal_t* tmp = this->geom[0].CellSize();
487 488 489 490 491 492 493 494
    
    Vector_t hr;
    for (int i = 0; i < 3; ++i)
        hr[i] = tmp[i];
    
    bunch_mp->setBaseLevelMeshSpacing(hr);
}

frey_m's avatar
frey_m committed
495 496 497 498 499 500

const Vector_t& AmrBoxLib::getMeshScaling() const {
    return meshScaling_m;
}


frey_m's avatar
frey_m committed
501
Vektor<int, 3> AmrBoxLib::getBaseLevelGridPoints() const {
frey_m's avatar
frey_m committed
502
    const Box_t& bx = this->geom[0].Domain();
503
    
frey_m's avatar
frey_m committed
504 505
    const AmrIntVect_t& low = bx.smallEnd();
    const AmrIntVect_t& high = bx.bigEnd();
506 507 508 509
    
    return Vektor<int, 3>(high[0] - low[0] + 1,
                          high[1] - low[1] + 1,
                          high[2] - low[2] + 1);
510 511 512
}


frey_m's avatar
frey_m committed
513 514
inline const int& AmrBoxLib::maxLevel() const {
    return this->max_level;
515 516 517
}


frey_m's avatar
frey_m committed
518 519
inline const int& AmrBoxLib::finestLevel() const {
    return this->finest_level;
520 521 522
}


523 524 525 526 527
inline double AmrBoxLib::getT() const {
    return bunch_mp->getT();
}


frey_m's avatar
frey_m committed
528
void AmrBoxLib::redistributeGrids(int /*how*/) {
529 530 531 532 533 534 535 536 537
//
//    // copied + modified version of AMReX_Amr.cpp
//    AmrProcMap_t::InitProximityMap();
////     AmrProcMap_t::Initialize();
//
//    AmrGridContainer_t allBoxes(finest_level + 1);
//    for(unsigned int ilev = 0; ilev < allBoxes.size(); ++ilev) {
//        allBoxes[ilev] = boxArray(ilev);
//    }
frey_m's avatar
frey_m committed
538
//    amrex::Vector<AmrIntArray_t> mLDM;
539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585
//        
//    switch ( how ) {
//        case RANK_ZERO:
//            mLDM = AmrProcMap_t::MultiLevelMapRandom(this->ref_ratio,
//                                                     allBoxes,
//                                                     maxGridSize(0),
//                                                     0/*maxRank*/,
//                                                     0/*minRank*/);
//            break;
//        case PFC:
//            mLDM = AmrProcMap_t::MultiLevelMapPFC(this->ref_ratio,
//                                                  allBoxes,
//                                                  maxGridSize(0));
//            break;
//        case RANDOM:
//            mLDM = AmrProcMap_t::MultiLevelMapRandom(this->ref_ratio,
//                                                     allBoxes,
//                                                     maxGridSize(0));
//            break;
//        case KNAPSACK:
//            mLDM = AmrProcMap_t::MultiLevelMapKnapSack(this->ref_ratio,
//                                                       allBoxes,
//                                                       maxGridSize(0));
//            break;
//        default:
//            *gmsg << "We didn't redistribute the grids." << endl;
//            return;
//    }
//        
//    for(unsigned int iMap = 0; iMap < mLDM.size(); ++iMap) {
//        AmrField_t::MoveAllFabs(mLDM[iMap]);
//    }
//    
//    /*
//     * particles need to know the BoxArray
//     * and DistributionMapping
//     */
//    for(unsigned int ilev = 0; ilev < allBoxes.size(); ++ilev) {
//        layout_mp->SetParticleBoxArray(ilev, this->grids[ilev]);
//        layout_mp->SetParticleDistributionMap(ilev, this->dmap[ilev]);
//    }
//    
//    for(unsigned int iMap = 0; iMap < mLDM.size(); ++iMap) {
//        rho_m[iMap]->MoveFabs(mLDM[iMap]);
//        phi_m[iMap]->MoveFabs(mLDM[iMap]);
//        efield_m[iMap]->MoveFabs(mLDM[iMap]);
//    }
frey_m's avatar
frey_m committed
586 587 588
}


frey_m's avatar
frey_m committed
589
void AmrBoxLib::RemakeLevel (int lev, AmrReal_t /*time*/,
frey_m's avatar
frey_m committed
590 591
                             const AmrGrid_t& new_grids,
                             const AmrProcMap_t& new_dmap)
592 593 594 595
{
    SetBoxArray(lev, new_grids);
    SetDistributionMap(lev, new_dmap);
    
frey_m's avatar
frey_m committed
596
    //                                                      #comp  #ghosts cells
597
    rho_m[lev].reset(new AmrField_t(new_grids, new_dmap,    1,     0));
frey_m's avatar
frey_m committed
598
    phi_m[lev].reset(new AmrField_t(new_grids, new_dmap,    1,     1));
599 600 601
    for (int j = 0; j < AMREX_SPACEDIM; ++j) {
        efield_m[lev][j].reset(new AmrField_t(new_grids, new_dmap, 1,     1));
    }
602 603
    
    // including nghost = 1
604
    rho_m[lev]->setVal(0.0, 0);
605
    phi_m[lev]->setVal(0.0, 1);
606 607 608 609
    
    for (int j = 0; j < AMREX_SPACEDIM; ++j) {
        efield_m[lev][j]->setVal(0.0, 1);
    }
frey_m's avatar
frey_m committed
610 611 612 613 614 615 616
    
    /*
     * particles need to know the BoxArray
     * and DistributionMapping
     */
     layout_mp->SetParticleBoxArray(lev, new_grids);
     layout_mp->SetParticleDistributionMap(lev, new_dmap);
frey_m's avatar
frey_m committed
617 618
     
     layout_mp->buildLevelMask(lev, this->nErrorBuf(lev));
619 620 621
}


frey_m's avatar
frey_m committed
622
void AmrBoxLib::MakeNewLevel (int lev, AmrReal_t /*time*/,
frey_m's avatar
frey_m committed
623 624
                              const AmrGrid_t& new_grids,
                              const AmrProcMap_t& new_dmap)
625 626 627 628
{
    SetBoxArray(lev, new_grids);
    SetDistributionMap(lev, new_dmap);
    
frey_m's avatar
frey_m committed
629
    //                                                      #comp  #ghosts cells
630
    rho_m[lev].reset(new AmrField_t(new_grids, new_dmap,    1,     0));
frey_m's avatar
frey_m committed
631
    phi_m[lev].reset(new AmrField_t(new_grids, new_dmap,    1,     1));
632 633 634 635
    
    for (int j = 0; j < AMREX_SPACEDIM; ++j) {
        efield_m[lev][j].reset(new AmrField_t(new_grids, new_dmap, 1,     1));
    }
636 637
    
    // including nghost = 1
638
    rho_m[lev]->setVal(0.0, 0);
639
    phi_m[lev]->setVal(0.0, 1);
640 641 642
    for (int j = 0; j < AMREX_SPACEDIM; ++j) {
        efield_m[lev][j]->setVal(0.0, 1);
    }
643 644
    
    
frey_m's avatar
frey_m committed
645 646 647 648 649 650 651
    
    /*
     * particles need to know the BoxArray
     * and DistributionMapping
     */
    layout_mp->SetParticleBoxArray(lev, new_grids);
    layout_mp->SetParticleDistributionMap(lev, new_dmap);
frey_m's avatar
frey_m committed
652 653
    
    layout_mp->buildLevelMask(lev, this->nErrorBuf(lev));
654 655 656 657
}


void AmrBoxLib::ClearLevel(int lev) {
frey_m's avatar
frey_m committed
658
    rho_m[lev].reset(nullptr);
659
    phi_m[lev].reset(nullptr);
660 661 662
    for (int j = 0; j < AMREX_SPACEDIM; ++j) {
        efield_m[lev][j].reset(nullptr);
    }
663 664
    ClearBoxArray(lev);
    ClearDistributionMap(lev);
665 666

    this->isFirstTagging_m[lev] = true;
667 668 669
}


frey_m's avatar
AMR:  
frey_m committed
670 671 672
void AmrBoxLib::ErrorEst(int lev, TagBoxArray_t& tags,
                         AmrReal_t time, int ngrow)
{
frey_m's avatar
frey_m committed
673
    *gmsg << level2 << "*         Start tagging of level " << lev << endl;
frey_m's avatar
frey_m committed
674
    
675
    switch ( tagging_m ) {
frey_m's avatar
AMR:  
frey_m committed
676
        case CHARGE_DENSITY:
677 678
            tagForChargeDensity_m(lev, tags, time, ngrow);
            break;
679
        case POTENTIAL:
680 681
            tagForPotentialStrength_m(lev, tags, time, ngrow);
            break;
682 683
        case EFIELD:
            tagForEfield_m(lev, tags, time, ngrow);
684
            break;
frey_m's avatar
frey_m committed
685 686 687 688 689 690 691 692 693
        case MOMENTA:
            tagForMomenta_m(lev, tags, time, ngrow);
            break;
        case MAX_NUM_PARTICLES:
            tagForMaxNumParticles_m(lev, tags, time, ngrow);
            break;
        case MIN_NUM_PARTICLES:
            tagForMinNumParticles_m(lev, tags, time, ngrow);
            break;
694 695 696 697
        default:
            tagForChargeDensity_m(lev, tags, time, ngrow);
            break;
    }
frey_m's avatar
frey_m committed
698
    
frey_m's avatar
frey_m committed
699
    *gmsg << level2 << "*         Finished tagging of level " << lev << endl;
700 701 702
}


frey_m's avatar
frey_m committed
703 704 705
void AmrBoxLib::MakeNewLevelFromScratch(int /*lev*/, AmrReal_t /*time*/,
                                  const AmrGrid_t& /*ba*/,
                                  const AmrProcMap_t& /*dm*/)
frey_m's avatar
AMR:  
frey_m committed
706 707 708 709 710
{
    throw OpalException("AmrBoxLib::MakeNewLevelFromScratch()", "Shouldn't be called.");
}


frey_m's avatar
frey_m committed
711 712 713
void AmrBoxLib::MakeNewLevelFromCoarse (int /*lev*/, AmrReal_t /*time*/,
                                        const AmrGrid_t& /*ba*/,
                                        const AmrProcMap_t& /*dm*/)
frey_m's avatar
AMR:  
frey_m committed
714 715 716 717 718
{
    throw OpalException("AmrBoxLib::MakeNewLevelFromCoarse()", "Shouldn't be called.");
}


719 720 721 722 723 724
void AmrBoxLib::doRegrid_m(int lbase, double time) {
    int new_finest = 0;
    AmrGridContainer_t new_grids(finest_level+2);
    
    MakeNewGrids(lbase, time, new_finest, new_grids);

snuverink_j's avatar
snuverink_j committed
725
    PAssert(new_finest <= finest_level+1);
726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754
    
    for (int lev = lbase+1; lev <= new_finest; ++lev)
    {
        if (lev <= finest_level) // an old level
        {
            if (new_grids[lev] != grids[lev]) // otherwise nothing
            {
                AmrProcMap_t new_dmap(new_grids[lev]);
                RemakeLevel(lev, time, new_grids[lev], new_dmap);
            }
        }
        else  // a new level
        {
            AmrProcMap_t new_dmap(new_grids[lev]);
            MakeNewLevel(lev, time, new_grids[lev], new_dmap);
        }
        
        layout_mp->setFinestLevel(new_finest);
    }

    // now we are safe to delete deprecated levels
    for (int lev = new_finest+1; lev <= finest_level; ++lev) {
        ClearLevel(lev);
    }
    
    finest_level = new_finest;
}


frey_m's avatar
frey_m committed
755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770
void AmrBoxLib::preRegrid_m() {
    /* In case of E-field or potential tagging
     * in combination of binning we need to make
     * sure we do not lose accuracy when tagging since
     * the grid data of the potential or e-field are only
     * non-zero for the last bin causing the tagging to refine
     * only there.
     * So, we need to solve the Poisson problem first assuming
     * a single bin only
     */
    if ( tagging_m == POTENTIAL || tagging_m == EFIELD ) {
        this->solvePoisson_m();
    }
}


771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809
void AmrBoxLib::postRegrid_m(int old_finest) {
    /* ATTENTION: In this call the bunch has to be updated
     * since each particle needs to be assigned to its latest
     * grid and sent to the corresponding MPI-process.
     * 
     * We need to update the bunch before we clear
     * levels, otherwise the particle level counter
     * is invalidated.
     */

    // redistribute particles and assign correct grid and level
    AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase();
    amrpbase_p->update(0, finest_level, true);

    // check if no particles left on higher levels
    const auto& LocalNumPerLevel = amrpbase_p->getLocalNumPerLevel();

    for (int lev = finest_level+1; lev <= old_finest; ++lev) {
        if ( LocalNumPerLevel.getLocalNumAtLevel(lev) != 0 ) {
            throw OpalException("AmrBoxLib::postRegrid_m()",
                                "Still particles on level " + std::to_string(lev) + "!");
        }
        layout_mp->ClearParticleBoxArray(lev);
        layout_mp->ClearParticleDistributionMap(lev);
        layout_mp->clearLevelMask(lev);
    }

    if ( bunch_mp->getLocalNum() != LocalNumPerLevel.getLocalNumUpToLevel(finest_level) ) {
        std::string localnum = std::to_string(bunch_mp->getLocalNum());
        std::string levelnum = std::to_string(LocalNumPerLevel.getLocalNumUpToLevel(finest_level));
        throw OpalException("AmrBoxLib::postRegrid_m()",
                            "Number of particles do not agree: " + localnum + " != " + levelnum);
    }

    PoissonSolver *solver = bunch_mp->getFieldSolver();
    solver->hasToRegrid();
}


frey_m's avatar
frey_m committed
810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858
double AmrBoxLib::solvePoisson_m() {
    AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase();

    /// from charge (C) to charge density (C/m^3).
    amrpbase_p->scatter(bunch_mp->Q, this->rho_m, bunch_mp->R, 0, finest_level, bunch_mp->Bin);

    int baseLevel = 0;

    // mesh scaling for solver
    meshScaling_m = Vector_t(1.0, 1.0, 1.0);

    // charge density is in rho_m
    // calculate Possion equation (with coefficient: -1/(eps))
    for (int i = 0; i <= finest_level; ++i) {
        if ( this->rho_m[i]->contains_nan(false) )
            throw OpalException("AmrBoxLib::solvePoisson_m() ",
                                "NANs at level " + std::to_string(i) + ".");
        this->rho_m[i]->mult(-1.0 / Physics::epsilon_0, 0, 1);
    }

    // find maximum and normalize each level (faster convergence)
    double l0norm = 1.0;
    for (int i = 0; i <= finest_level; ++i)
        l0norm = std::max(l0norm, this->rho_m[i]->norm0(0));

    for (int i = 0; i <= finest_level; ++i) {
        this->rho_m[i]->mult(1.0 / l0norm, 0, 1);
    }

    PoissonSolver *solver = bunch_mp->getFieldSolver();

    IpplTimings::startTimer(this->amrSolveTimer_m);
    solver->solve(rho_m, phi_m, efield_m, baseLevel, finest_level);
    IpplTimings::stopTimer(this->amrSolveTimer_m);

    // make sure ghost cells are filled
    for (int i = 0; i <= finest_level; ++i) {
        phi_m[i]->FillBoundary(geom[i].periodicity());
        for (int j = 0; j < AMREX_SPACEDIM; ++j) {
            efield_m[i][j]->FillBoundary(geom[i].periodicity());
        }
    }

    this->fillPhysbc_m(*(this->phi_m[0]), 0);

    return l0norm;
}


frey_m's avatar
AMR:  
frey_m committed
859
void AmrBoxLib::tagForChargeDensity_m(int lev, TagBoxArray_t& tags,
frey_m's avatar
frey_m committed
860
                                      AmrReal_t /*time*/, int /*ngrow*/)
frey_m's avatar
AMR:  
frey_m committed
861
{
862
    // we need to update first
frey_m's avatar
AMR:  
frey_m committed
863
    AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase();
864 865
    amrpbase_p->update(0, lev, true);
    
frey_m's avatar
frey_m committed
866
    for (int i = lev; i <= finest_level; ++i)
frey_m's avatar
frey_m committed
867
        rho_m[i]->setVal(0.0, rho_m[i]->nGrow());
frey_m's avatar
frey_m committed
868
    
frey_m's avatar
frey_m committed
869
    // the new scatter function averages the value also down to the coarsest level
870
    amrpbase_p->scatter(bunch_mp->Q, rho_m,
871
                        bunch_mp->R, 0, lev, bunch_mp->Bin);
frey_m's avatar
frey_m committed
872
    
frey_m's avatar
frey_m committed
873
    const double& scalefactor = amrpbase_p->getScalingFactor();
874 875 876 877
    
    for (int i = lev; i <= finest_level; ++i)
        rho_m[i]->mult(scalefactor, 0, 1);
    
frey_m's avatar
frey_m committed
878 879
    const int clearval = TagBox_t::CLEAR;
    const int   tagval = TagBox_t::SET;
frey_m's avatar
frey_m committed
880 881 882 883 884

#ifdef _OPENMP
#pragma omp parallel
#endif
    {
frey_m's avatar
frey_m committed
885 886
        for (MFIter_t mfi(*rho_m[lev], true); mfi.isValid(); ++mfi) {
            const Box_t&  tilebx  = mfi.tilebox();
frey_m's avatar
frey_m committed
887
            
frey_m's avatar
frey_m committed
888
            TagBox_t&     tagfab  = tags[mfi];
frey_m's avatar
frey_m committed
889
            FArrayBox_t&  fab = (*rho_m[lev])[mfi];
frey_m's avatar
frey_m committed
890 891 892
            
            const int*  tlo     = tilebx.loVect();
            const int*  thi     = tilebx.hiVect();
frey_m's avatar
frey_m committed
893 894 895 896 897 898 899 900 901 902 903 904 905 906
            
            for (int i = tlo[0]; i <= thi[0]; ++i) {
                for (int j = tlo[1]; j <= thi[1]; ++j) {
                    for (int k = tlo[2]; k <= thi[2]; ++k) {
                        
                        amrex::IntVect iv(D_DECL(i,j,k));
                        
                        if ( std::abs( fab(iv) ) >= chargedensity_m )
                            tagfab(iv) = tagval;
                        else
                            tagfab(iv) = clearval;
                    }
                }
            }
frey_m's avatar
frey_m committed
907 908
        }
    }
909 910 911
}


frey_m's avatar
frey_m committed
912 913 914
void AmrBoxLib::tagForPotentialStrength_m(int lev, TagBoxArray_t& tags,
                                          AmrReal_t time, int ngrow)
{
frey_m's avatar
frey_m committed
915
    /* Tag all cells for refinement
frey_m's avatar
frey_m committed
916 917 918
     * where the value of the potential is higher than 75 percent of the maximum potential
     * value of this level.
     */
frey_m's avatar
frey_m committed
919
    if ( !isPoissonSolved_m || !phi_m[lev]->ok() || this->isFirstTagging_m[lev] ) {
frey_m's avatar
frey_m committed
920
        *gmsg << level2 << "* Level " << lev << ": We need to perform "
921
              << "charge tagging if a new level is created." << endl;
922
        this->tagForChargeDensity_m(lev, tags, time, ngrow);
frey_m's avatar
frey_m committed
923
        this->isFirstTagging_m[lev] = false;
924
        return;
frey_m's avatar
frey_m committed
925
    }
frey_m's avatar
frey_m committed
926
    
927
    // tag cells for refinement
frey_m's avatar
frey_m committed
928 929
    const int clearval = TagBox_t::CLEAR;
    const int   tagval = TagBox_t::SET;
frey_m's avatar
frey_m committed
930

frey_m's avatar
frey_m committed
931
    AmrReal_t threshold = phi_m[lev]->norm0(0, 0 /*nghost*/, false /*local*/);
frey_m's avatar
frey_m committed
932
    
933
    threshold *= scaling_m;
frey_m's avatar
frey_m committed
934 935 936 937 938
    
#ifdef _OPENMP
#pragma omp parallel
#endif
    {
frey_m's avatar
frey_m committed
939
        for (MFIter_t mfi(*phi_m[lev], true); mfi.isValid(); ++mfi) {
frey_m's avatar
frey_m committed
940
            
frey_m's avatar
frey_m committed
941
            const Box_t&  tilebx  = mfi.tilebox();
frey_m's avatar
frey_m committed
942
            TagBox_t&     tagfab  = tags[mfi];
frey_m's avatar
frey_m committed
943
            FArrayBox_t&  fab     = (*phi_m[lev])[mfi];
frey_m's avatar
frey_m committed
944 945 946
            
            const int*  tlo     = tilebx.loVect();
            const int*  thi     = tilebx.hiVect();
frey_m's avatar
frey_m committed
947 948 949 950 951 952 953 954 955 956 957 958 959 960
            
            for (int i = tlo[0]; i <= thi[0]; ++i) {
                for (int j = tlo[1]; j <= thi[1]; ++j) {
                    for (int k = tlo[2]; k <= thi[2]; ++k) {
                        
                        amrex::IntVect iv(D_DECL(i,j,k));
                        
                        if ( std::abs( fab(iv) ) >= threshold )
                            tagfab(iv) = tagval;
                        else
                            tagfab(iv) = clearval;
                    }
                }
            }
frey_m's avatar
frey_m committed
961 962
        }
    }
963 964 965
}


frey_m's avatar
frey_m committed
966 967 968
void AmrBoxLib::tagForEfield_m(int lev, TagBoxArray_t& tags,
                               AmrReal_t time, int ngrow)
{
frey_m's avatar
frey_m committed
969
    /* Tag all cells for refinement
970
     * where the value of the efield is higher than 75 percent of the maximum efield
frey_m's avatar
frey_m committed
971 972
     * value of this level.
     */
frey_m's avatar
frey_m committed
973
    if ( !isPoissonSolved_m || !efield_m[lev][0]->ok() || this->isFirstTagging_m[lev] ) {
frey_m's avatar
frey_m committed
974
        *gmsg << level2 << "* Level " << lev << ": We need to perform "
975
              << "charge tagging if a new level is created." << endl;
976
        this->tagForChargeDensity_m(lev, tags, time, ngrow);
frey_m's avatar
frey_m committed
977
        this->isFirstTagging_m[lev] = false;
frey_m's avatar
frey_m committed
978
        return;
frey_m's avatar
frey_m committed
979
    }
frey_m's avatar
frey_m committed
980
    
981
    // tag cells for refinement
frey_m's avatar
frey_m committed
982 983
    const int clearval = TagBox_t::CLEAR;
    const int   tagval = TagBox_t::SET;
frey_m's avatar
frey_m committed
984
    
frey_m's avatar
frey_m committed
985
    // obtain maximum absolute value
986 987 988 989 990 991
    amrex::Vector<AmrReal_t> threshold(AMREX_SPACEDIM);
    
    for (int i = 0; i < 3; ++i) {
        threshold[i] = efield_m[lev][i]->norm0(0,
                                            0 /*nghosts*/,
                                            false /*local*/);
frey_m's avatar
frey_m committed
992 993
    
        threshold[i] *= scaling_m;
994
    }
frey_m's avatar
frey_m committed
995 996 997 998 999
    
#ifdef _OPENMP
#pragma omp parallel
#endif
    {
1000
        for (MFIter_t mfi(this->grids[lev], this->dmap[lev], true); mfi.isValid(); ++mfi) {
frey_m's avatar
frey_m committed
1001
            
frey_m's avatar
frey_m committed
1002
            const Box_t&  tilebx  = mfi.tilebox();
frey_m's avatar
frey_m committed
1003
            TagBox_t&     tagfab  = tags[mfi];
1004 1005 1006
            FArrayBox_t&  xfab     = (*efield_m[lev][0])[mfi];
            FArrayBox_t&  yfab     = (*efield_m[lev][1])[mfi];
            FArrayBox_t&  zfab     = (*efield_m[lev][2])[mfi];
frey_m's avatar
frey_m committed
1007

frey_m's avatar
frey_m committed
1008 1009 1010 1011 1012 1013
            const int*  tlo     = tilebx.loVect();
            const int*  thi     = tilebx.hiVect();

            for (int i = tlo[0]; i <= thi[0]; ++i) {
                for (int j = tlo[1]; j <= thi[1]; ++j) {
                    for (int k = tlo[2]; k <= thi[2]; ++k) {
frey_m's avatar
frey_m committed
1014
                        AmrIntVect_t iv(i,j,k);
1015
                        if (std::abs(xfab(iv)) >= threshold[0])
frey_m's avatar
frey_m committed
1016
                            tagfab(iv) = tagval;
1017
                        else if (std::abs(yfab(iv)) >= threshold[1])
frey_m's avatar
frey_m committed
1018
                            tagfab(iv) = tagval;
1019
                        else if (std::abs(zfab(iv)) >= threshold[2])
frey_m's avatar
frey_m committed
1020 1021 1022 1023 1024 1025
                            tagfab(iv) = tagval;
                        else
                            tagfab(iv) = clearval;
                    }
                }
            }
frey_m's avatar
frey_m committed
1026 1027 1028 1029 1030 1031
        }
    }
}


void AmrBoxLib::tagForMomenta_m(int lev, TagBoxArray_t& tags,
frey_m's avatar
frey_m committed
1032
                                AmrReal_t /*time*/, int /*ngrow*/)
frey_m's avatar
frey_m committed
1033 1034
{
    AmrPartBunch::pbase_t* amrpbase_p = bunch_mp->getAmrParticleBase();