AmrMultiGrid.cpp 75 KB
Newer Older
1 2 3 4
//
// Class AmrMultiGrid
//   Main class of the AMR Poisson multigrid solver.
//   It implements the multigrid solver described in https://doi.org/10.1016/j.cpc.2019.106912
5 6 7 8
//
// Copyright (c) 2017 - 2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
9 10 11 12 13 14 15 16 17 18 19 20 21
// 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/>.
//
22 23
#include "AmrMultiGrid.h"

frey_m's avatar
frey_m committed
24
#include <algorithm>
frey_m's avatar
frey_m committed
25
#include <functional>
26
#include <map>
frey_m's avatar
frey_m committed
27
#include <numeric>
28

29 30
#include "OPALconfig.h"
#include "AbstractObjects/OpalData.h"
31
#include "Utilities/OpalException.h"
32
#include "Utilities/Timer.h"
33
#include "Utilities/Util.h"
frey_m's avatar
frey_m committed
34
#include <AMReX_ParallelDescriptor.H>
35

36 37
#include <boost/filesystem.hpp>

38 39 40 41
#if AMR_MG_WRITE
    #include <iomanip>
#endif

42 43
AmrMultiGrid::AmrMultiGrid(AmrBoxLib* itsAmrObject_p,
                           const std::string& bsolver,
44
                           const std::string& prec,
45
                           const bool& rebalance,
frey_m's avatar
frey_m committed
46
                           const std::string& reuse,
47 48 49 50 51 52 53
                           const std::string& bcx,
                           const std::string& bcy,
                           const std::string& bcz,
                           const std::string& smoother,
                           const std::size_t& nSweeps,
                           const std::string& interp,
                           const std::string& norm)
frey_m's avatar
frey_m committed
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
    : AmrPoissonSolver<AmrBoxLib>(itsAmrObject_p)
    , comm_mp( new comm_t( amrex::ParallelDescriptor::Communicator() ) )
    , nIter_m(0)
    , bIter_m(0)
    , maxiter_m(100)
    , nSweeps_m(nSweeps)
    , mglevel_m(0)
    , lbase_m(0)
    , lfine_m(0)
    , nlevel_m(1)
    , nBcPoints_m(0)
    , eps_m(1.0e-10)
    , verbose_m(false)
    , fname_m(OpalData::getInstance()->getInputBasename() + std::string(".solver"))
    , flag_m(std::ios::out)
69 70
{
    node_mp = KokkosClassic::Details::getNode<amr::node_t>(); //KokkosClassic::DefaultNode::getDefaultNode();
kraus's avatar
kraus committed
71

72
#if AMR_MG_TIMER
73
    this->initTimer_m();
74
#endif
kraus's avatar
kraus committed
75

76 77 78 79 80
    const Boundary bcs[AMREX_SPACEDIM] = {
        this->convertToEnumBoundary_m(bcx),
        this->convertToEnumBoundary_m(bcy),
        this->convertToEnumBoundary_m(bcz)
    };
kraus's avatar
kraus committed
81

82
    this->initPhysicalBoundary_m(&bcs[0]);
kraus's avatar
kraus committed
83

84
    smootherType_m = this->convertToEnumSmoother_m(smoother);
kraus's avatar
kraus committed
85

86
    norm_m = this->convertToEnumNorm_m(norm);
kraus's avatar
kraus committed
87

88 89
    const Interpolater interpolater = this->convertToEnumInterpolater_m(interp);
    this->initInterpolater_m(interpolater);
kraus's avatar
kraus committed
90

91
    // interpolater for crse-fine-interface
92
    this->initCrseFineInterp_m(Interpolater::LAGRANGE);
kraus's avatar
kraus committed
93

94
    // preconditioner
95
    const Preconditioner precond = this->convertToEnumPreconditioner_m(prec);
frey_m's avatar
frey_m committed
96
    this->initPrec_m(precond, rebalance, reuse);
kraus's avatar
kraus committed
97

98
    // base level solver
99
    const BaseSolver solver = this->convertToEnumBaseSolver_m(bsolver);
frey_m's avatar
frey_m committed
100
    this->initBaseSolver_m(solver, rebalance, reuse);
kraus's avatar
kraus committed
101

102 103 104 105 106 107
    if (boost::filesystem::exists(fname_m)) {
        flag_m = std::ios::app;
        INFOMSG("Appending solver information to existing file: " << fname_m << endl);
    } else {
        INFOMSG("Creating new file for solver information: " << fname_m << endl);
    }
108 109
}

110

111 112 113
void AmrMultiGrid::solve(AmrScalarFieldContainer_t &rho,
                         AmrScalarFieldContainer_t &phi,
                         AmrVectorFieldContainer_t &efield,
114 115 116 117 118 119 120
                         unsigned short baseLevel,
                         unsigned short finestLevel,
                         bool prevAsGuess)
{
    lbase_m = baseLevel;
    lfine_m = finestLevel;
    nlevel_m = lfine_m - lbase_m + 1;
kraus's avatar
kraus committed
121

frey_m's avatar
frey_m committed
122 123
    /* we cannot use the previous solution
     * if we have to regrid (AmrPoissonSolver::hasToRegrid())
kraus's avatar
kraus committed
124
     *
frey_m's avatar
frey_m committed
125 126
     * regrid_m is set in AmrBoxlib::regrid()
     */
frey_m's avatar
frey_m committed
127 128 129 130
    bool reset = !prevAsGuess;

    if ( this->regrid_m )
        reset = true;
frey_m's avatar
frey_m committed
131

frey_m's avatar
frey_m committed
132
    this->initLevels_m(rho, itsAmrObject_mp->Geom(), this->regrid_m);
frey_m's avatar
frey_m committed
133

134
    // build all necessary matrices and vectors
frey_m's avatar
frey_m committed
135
    this->setup_m(rho, phi, this->regrid_m);
frey_m's avatar
frey_m committed
136

frey_m's avatar
frey_m committed
137
    this->initGuess_m(reset);
frey_m's avatar
frey_m committed
138

139
    // actual solve
140
    scalar_t error = this->iterate_m();
frey_m's avatar
frey_m committed
141

frey_m's avatar
frey_m committed
142 143 144
    for (int lev = nlevel_m - 1; lev > -1; --lev) {
        averageDown_m(lev);
    }
kraus's avatar
kraus committed
145

146
    // write efield to AMReX
kraus's avatar
kraus committed
147
    this->computeEfield_m(efield);
148
    
149 150 151
    // copy solution back
    for (int lev = 0; lev < nlevel_m; ++lev) {
        int ilev = lbase_m + lev;
kraus's avatar
kraus committed
152

frey_m's avatar
frey_m committed
153
        phi[ilev]->setVal(0.0, phi[ilev]->nGrow());
kraus's avatar
kraus committed
154

155
        this->trilinos2amrex_m(lev, 0, *phi[ilev], mglevel_m[lev]->phi_p);
156
    }
kraus's avatar
kraus committed
157

158 159
    if ( verbose_m )
        this->writeSDDSData_m(error);
kraus's avatar
kraus committed
160

161 162
    // we can now reset
    this->regrid_m = false;
163 164 165
}


166 167 168 169
void AmrMultiGrid::setNumberOfSweeps(const std::size_t& nSweeps) {
    if ( nSweeps < 0 )
        throw OpalException("AmrMultiGrid::setNumberOfSweeps()",
                            "The number of smoothing sweeps needs to be non-negative!");
kraus's avatar
kraus committed
170

171 172 173 174 175 176 177 178
    nSweeps_m = nSweeps;
}


void AmrMultiGrid::setMaxNumberOfIterations(const std::size_t& maxiter) {
    if ( maxiter < 1 )
        throw OpalException("AmrMultiGrid::setMaxNumberOfIterations()",
                            "The max. number of iterations needs to be positive!");
kraus's avatar
kraus committed
179

180 181 182 183
    maxiter_m = maxiter;
}


184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
std::size_t AmrMultiGrid::getNumIters() {
    return nIter_m;
}


AmrMultiGrid::scalar_t AmrMultiGrid::getLevelResidualNorm(lo_t level) {
    return evalNorm_m(mglevel_m[level]->residual_p);
}


void AmrMultiGrid::setVerbose(bool verbose) {
    verbose_m = verbose;
}


199 200 201 202
void AmrMultiGrid::initPhysicalBoundary_m(const Boundary* bc)
{
    // make sure it's reset
    nBcPoints_m = 0;
kraus's avatar
kraus committed
203

204 205
    for (unsigned int i = 0; i < AMREX_SPACEDIM; ++i) {
        switch ( bc[i] ) {
206
            case Boundary::DIRICHLET:
207
                bc_m[i].reset( new AmrDirichletBoundary<AmrMultiGridLevel_t>() );
208 209
                break;
            case Boundary::OPEN:
210
                bc_m[i].reset( new AmrOpenBoundary<AmrMultiGridLevel_t>() );
211 212
                break;
            case Boundary::PERIODIC:
213
                bc_m[i].reset( new AmrPeriodicBoundary<AmrMultiGridLevel_t>() );
214 215
                break;
            default:
216 217 218 219
                throw OpalException("AmrMultiGrid::initPhysicalBoundary_m()",
                                    "This type of boundary is not supported");
        }
        // we use the maximum in order to build matrices
frey_m's avatar
frey_m committed
220
        go_t tmp = bc_m[i]->getNumberOfPoints();
221 222
        if ( nBcPoints_m < tmp )
            nBcPoints_m = tmp;
223 224 225 226
    }
}


frey_m's avatar
frey_m committed
227 228
void AmrMultiGrid::initLevels_m(const amrex::Vector<AmrField_u>& rho,
                                const amrex::Vector<AmrGeometry_t>& geom,
frey_m's avatar
frey_m committed
229
                                bool regrid)
230
{
frey_m's avatar
frey_m committed
231
    if ( !regrid )
232
        return;
kraus's avatar
kraus committed
233

234 235 236 237 238
    // although we do a resize afterwards, we do this to be safe
    for (int lev = nlevel_m; lev < (int)mglevel_m.size(); ++lev) {
        mglevel_m[lev].reset(nullptr);
    }

239
    mglevel_m.resize(nlevel_m);
kraus's avatar
kraus committed
240

frey_m's avatar
frey_m committed
241
    amrex::Periodicity period(AmrIntVect_t(D_DECL(0, 0, 0)));
kraus's avatar
kraus committed
242

243
    AmrIntVect_t rr = AmrIntVect_t(D_DECL(2, 2, 2));
kraus's avatar
kraus committed
244

245
    for (int lev = 0; lev < nlevel_m; ++lev) {
246
        int ilev = lbase_m + lev;
kraus's avatar
kraus committed
247

frey_m's avatar
frey_m committed
248 249 250 251 252 253 254 255 256 257 258 259 260
        // do not initialize base level every time
        if (mglevel_m[lev] == nullptr || lev > lbase_m) {
            mglevel_m[lev].reset(new AmrMultiGridLevel_t(itsAmrObject_mp->getMeshScaling(),
                                                         rho[ilev]->boxArray(),
                                                         rho[ilev]->DistributionMap(),
                                                         geom[ilev],
                                                         rr,
                                                         bc_m,
                                                         comm_mp,
                                                         node_mp));
        } else {
            mglevel_m[lev]->buildLevelMask();
        }
261 262 263 264 265 266

        mglevel_m[lev]->refmask.reset(
            new AmrMultiGridLevel_t::mask_t(mglevel_m[lev]->grids,
                                            mglevel_m[lev]->dmap, 1, 2)
            );
        mglevel_m[lev]->refmask->setVal(AmrMultiGridLevel_t::Refined::NO, 2);
frey_m's avatar
frey_m committed
267
        mglevel_m[lev]->refmask->FillBoundary(period);
268 269 270 271 272 273 274

        amrex::BoxArray ba = mglevel_m[lev]->grids;
        ba.coarsen(rr);
        mglevel_m[lev]->crsemask.reset(
            new AmrMultiGridLevel_t::mask_t(ba,
                                            mglevel_m[lev]->dmap, 1, 2)
            );
275
        mglevel_m[lev]->crsemask->setVal(AmrMultiGridLevel_t::Refined::NO, 2);
276 277 278 279 280 281 282 283 284
    }

    for (int lev = 1; lev < nlevel_m; ++lev) {
        mglevel_m[lev]->crsemask->setVal(AmrMultiGridLevel_t::Refined::YES, 0);

        // used for boundary interpolation --> replaces expensive calls to isBoundary
        mglevel_m[lev]->crsemask->setDomainBndry(AmrMultiGridLevel_t::Mask::PHYSBNDRY,
                                                 mglevel_m[lev-1]->geom); //FIXME: geometry of lev - 1
        // really needed ?
frey_m's avatar
frey_m committed
285
        mglevel_m[lev]->crsemask->FillBoundary(period);
286 287 288 289 290 291 292 293 294 295 296 297 298
    }

    /* to complete initialization we need to fill
     * the mask of refinement
     */
    for (int lev = 0; lev < nlevel_m-1; ++lev) {
        // get boxarray with refined cells
        amrex::BoxArray ba = mglevel_m[lev]->grids;
        ba.refine(rr);
        ba = amrex::intersect(mglevel_m[lev+1]->grids, ba);
        ba.coarsen(rr);

        // refined cells
frey_m's avatar
frey_m committed
299
        amrex::DistributionMapping dmap(ba, comm_mp->getSize());
300 301 302 303 304 305 306 307 308 309 310 311
        AmrMultiGridLevel_t::mask_t refined(ba, dmap, 1, 0);
        refined.setVal(AmrMultiGridLevel_t::Refined::YES);
//      refined.setDomainBndry(AmrMultiGridLevel_t::Mask::PHYSBNDRY, mglevel_m[lev]->geom);

        // fill intersection with YES
        mglevel_m[lev]->refmask->copy(refined, 0, 0, 1, 0, 2);

        /* physical boundary cells will never be refined cells
         * since they are ghost cells
         */
        mglevel_m[lev]->refmask->setDomainBndry(AmrMultiGridLevel_t::Mask::PHYSBNDRY,
                                                mglevel_m[lev]->geom);
kraus's avatar
kraus committed
312

frey_m's avatar
frey_m committed
313
        mglevel_m[lev]->refmask->FillBoundary(period);
314 315 316 317 318 319 320 321 322
    }
}


void AmrMultiGrid::clearMasks_m() {
    for (int lev = 0; lev < nlevel_m; ++lev) {
        mglevel_m[lev]->refmask.reset(nullptr);
        mglevel_m[lev]->crsemask.reset(nullptr);
        mglevel_m[lev]->mask.reset(nullptr);
323 324 325 326
    }
}


frey_m's avatar
frey_m committed
327 328
void AmrMultiGrid::initGuess_m(bool reset) {
    if ( !reset )
329
        return;
kraus's avatar
kraus committed
330

331 332 333
    // reset
    for (int lev = 0; lev < nlevel_m; ++lev)
        mglevel_m[lev]->phi_p->putScalar(0.0);
334 335 336
}


337
AmrMultiGrid::scalar_t AmrMultiGrid::iterate_m() {
kraus's avatar
kraus committed
338

339
    // initial error
frey_m's avatar
frey_m committed
340 341
    std::vector<scalar_t> rhsNorms;
    std::vector<scalar_t> resNorms;
kraus's avatar
kraus committed
342

frey_m's avatar
frey_m committed
343
    this->initResidual_m(rhsNorms, resNorms);
kraus's avatar
kraus committed
344

frey_m's avatar
frey_m committed
345
    std::for_each(rhsNorms.begin(), rhsNorms.end(),
346
                  [this](double& val){ val *= eps_m; });
kraus's avatar
kraus committed
347

348
    nIter_m = 0;
349
    bIter_m = 0;
kraus's avatar
kraus committed
350

frey_m's avatar
frey_m committed
351
    while ( !isConverged_m(rhsNorms, resNorms) && nIter_m < maxiter_m ) {
kraus's avatar
kraus committed
352

353
        relax_m(lfine_m);
kraus's avatar
kraus committed
354

frey_m's avatar
frey_m committed
355 356 357 358 359 360 361 362
//         /* in contrast to algorithm, we average down now
//          * --> potential is valid also on coarse covered
//          * cells
//          * --> however, it may take 1-2 iterations longer
//          */
//         for (int lev = nlevel_m - 1; lev > -1; --lev) {
//             averageDown_m(lev);
//         }
kraus's avatar
kraus committed
363

364 365
        // update residual
        for (int lev = 0; lev < nlevel_m; ++lev) {
kraus's avatar
kraus committed
366

367 368
            this->residual_m(lev,
                             mglevel_m[lev]->residual_p,
369
                             mglevel_m[lev]->rho_p,
370
                             mglevel_m[lev]->phi_p);
371
        }
kraus's avatar
kraus committed
372 373


frey_m's avatar
frey_m committed
374 375
        for (lo_t lev = 0; lev < nlevel_m; ++lev)
            resNorms[lev] = getLevelResidualNorm(lev);
kraus's avatar
kraus committed
376

377
        ++nIter_m;
kraus's avatar
kraus committed
378

frey_m's avatar
frey_m committed
379
#if AMR_MG_WRITE
frey_m's avatar
frey_m committed
380 381
        this->writeResidualNorm_m();
#endif
kraus's avatar
kraus committed
382

383
        bIter_m += solver_mp->getNumIters();
384
    }
kraus's avatar
kraus committed
385

frey_m's avatar
frey_m committed
386 387 388
    return std::accumulate(resNorms.begin(),
                           resNorms.end(), 0.0,
                           std::plus<scalar_t>());
389 390 391
}


frey_m's avatar
frey_m committed
392 393 394 395
bool AmrMultiGrid::isConverged_m(std::vector<scalar_t>& rhsNorms,
                                 std::vector<scalar_t>& resNorms)
{
    return std::equal(resNorms.begin(), resNorms.end(),
frey_m's avatar
frey_m committed
396
                      rhsNorms.begin(),
frey_m's avatar
frey_m committed
397
                      std::less<scalar_t>());
frey_m's avatar
frey_m committed
398 399 400
}


401 402
void AmrMultiGrid::residual_m(const lo_t& level,
                              Teuchos::RCP<vector_t>& r,
403
                              const Teuchos::RCP<vector_t>& b,
404
                              const Teuchos::RCP<vector_t>& x)
405 406 407 408 409
{
    /*
     * r = b - A*x
     */
    if ( level < lfine_m ) {
kraus's avatar
kraus committed
410

411
        vector_t fine2crse(mglevel_m[level]->Awf_p->getDomainMap(), true);
kraus's avatar
kraus committed
412 413

        // get boundary for
414 415 416
        if ( mglevel_m[level]->Bfine_p != Teuchos::null ) {
            mglevel_m[level]->Bfine_p->apply(*mglevel_m[level+1]->phi_p, fine2crse);
        }
kraus's avatar
kraus committed
417

418 419 420 421 422
        // operation: fine2crse += A * x
        mglevel_m[level]->Awf_p->apply(*x, fine2crse,
                                       Teuchos::NO_TRANS,
                                       scalar_t(1.0),
                                       scalar_t(1.0));
kraus's avatar
kraus committed
423

424
        if ( mglevel_m[level]->Bcrse_p != Teuchos::null ) {
425 426 427 428 429 430
            // operation: fine2crse += B * phi^(l-1)
            mglevel_m[level]->Bcrse_p->apply(*mglevel_m[level-1]->phi_p,
                                             fine2crse,
                                             Teuchos::NO_TRANS,
                                             scalar_t(1.0),
                                             scalar_t(1.0));
431
        }
kraus's avatar
kraus committed
432

frey_m's avatar
frey_m committed
433 434
        Teuchos::RCP<vector_t> uncov_Ax = Teuchos::rcp( new vector_t(mglevel_m[level]->map_p, true) );
        mglevel_m[level]->UnCovered_p->apply(fine2crse, *uncov_Ax);
kraus's avatar
kraus committed
435

frey_m's avatar
frey_m committed
436
        Teuchos::RCP<vector_t> uncov_b = Teuchos::rcp( new vector_t(mglevel_m[level]->map_p, true) );
kraus's avatar
kraus committed
437

frey_m's avatar
frey_m committed
438
        mglevel_m[level]->UnCovered_p->apply(*b, *uncov_b);
kraus's avatar
kraus committed
439

440 441
        // ONLY subtract coarse rho
//         mglevel_m[level]->residual_p->putScalar(0.0);
kraus's avatar
kraus committed
442

frey_m's avatar
frey_m committed
443
        r->update(1.0, *uncov_b, -1.0, *uncov_Ax, 0.0);
kraus's avatar
kraus committed
444

445 446 447
    } else {
        /* finest level: Awf_p == Anf_p
         */
frey_m's avatar
frey_m committed
448
        Teuchos::RCP<vector_t> Ax = Teuchos::rcp( new vector_t(mglevel_m[level]->map_p, true) );
449 450
        mglevel_m[level]->Anf_p->apply(*x, *Ax);
        
451
        if ( mglevel_m[level]->Bcrse_p != Teuchos::null ) {
frey_m's avatar
frey_m committed
452
            // operationr: Ax += B * phi^(l-1)
453
            mglevel_m[level]->Bcrse_p->apply(*mglevel_m[level-1]->phi_p,
frey_m's avatar
frey_m committed
454
                                             *Ax,
455 456 457
                                             Teuchos::NO_TRANS,
                                             scalar_t(1.0),
                                             scalar_t(1.0));
458
        }
frey_m's avatar
frey_m committed
459
        r->update(1.0, *b, -1.0, *Ax, 0.0);
460 461 462 463
    }
}


464
void AmrMultiGrid::relax_m(const lo_t& level) {
kraus's avatar
kraus committed
465

466
    if ( level == lfine_m ) {
kraus's avatar
kraus committed
467

468 469 470
        if ( level == lbase_m ) {
            /* Anf_p == Awf_p
             */
471 472 473 474
            Teuchos::RCP<vector_t> Ax = Teuchos::rcp( new vector_t(mglevel_m[level]->map_p, true) );
            mglevel_m[level]->Anf_p->apply(*mglevel_m[level]->phi_p, *Ax);
            mglevel_m[level]->residual_p->update(1.0, *mglevel_m[level]->rho_p, -1.0, *Ax, 0.0);
            
475
        } else {
476 477
            this->residual_no_fine_m(level,
                                     mglevel_m[level]->residual_p,
478 479
                                     mglevel_m[level]->phi_p,
                                     mglevel_m[level-1]->phi_p,
480
                                     mglevel_m[level]->rho_p);
481 482
        }
    }
kraus's avatar
kraus committed
483

484
    if ( level > 0 ) {
kraus's avatar
kraus committed
485
        // phi^(l, save) = phi^(l)
486 487
        Teuchos::RCP<vector_t> phi_save = Teuchos::rcp( new vector_t(mglevel_m[level]->map_p) );
        Tpetra::deep_copy(*phi_save, *mglevel_m[level]->phi_p);
kraus's avatar
kraus committed
488

489
        mglevel_m[level-1]->error_p->putScalar(0.0);
kraus's avatar
kraus committed
490

491
        // smoothing
492 493 494
        this->smooth_m(level,
                       mglevel_m[level]->error_p,
                       mglevel_m[level]->residual_p);
kraus's avatar
kraus committed
495 496


497 498
        // phi = phi + e
        mglevel_m[level]->phi_p->update(1.0, *mglevel_m[level]->error_p, 1.0);
kraus's avatar
kraus committed
499

500 501 502 503
        /*
         * restrict
         */
        this->restrict_m(level);
kraus's avatar
kraus committed
504

505
        this->relax_m(level - 1);
kraus's avatar
kraus committed
506

507 508 509 510
        /*
         * prolongate / interpolate
         */
        this->prolongate_m(level);
kraus's avatar
kraus committed
511

512 513
        // residual update
        Teuchos::RCP<vector_t> tmp = Teuchos::rcp( new vector_t(mglevel_m[level]->map_p) );
514
        this->residual_no_fine_m(level, tmp,
515 516
                                 mglevel_m[level]->error_p,
                                 mglevel_m[level-1]->error_p,
517
                                 mglevel_m[level]->residual_p);
kraus's avatar
kraus committed
518

519
        Tpetra::deep_copy(*mglevel_m[level]->residual_p, *tmp);
kraus's avatar
kraus committed
520

521
        // delta error
frey_m's avatar
frey_m committed
522
        Teuchos::RCP<vector_t> derror = Teuchos::rcp( new vector_t(mglevel_m[level]->map_p, true) );
kraus's avatar
kraus committed
523

524
        // smoothing
525
        this->smooth_m(level, derror, mglevel_m[level]->residual_p);
kraus's avatar
kraus committed
526

527 528
        // e^(l) += de^(l)
        mglevel_m[level]->error_p->update(1.0, *derror, 1.0);
kraus's avatar
kraus committed
529

530 531
        // phi^(l) = phi^(l, save) + e^(l)
        mglevel_m[level]->phi_p->update(1.0, *phi_save, 1.0, *mglevel_m[level]->error_p, 0.0);
kraus's avatar
kraus committed
532

533 534 535 536 537
    } else {
        // e = A^(-1)r
#if AMR_MG_TIMER
    IpplTimings::startTimer(bottomTimer_m);
#endif
kraus's avatar
kraus committed
538

539 540
        solver_mp->solve(mglevel_m[level]->error_p,
                         mglevel_m[level]->residual_p);
kraus's avatar
kraus committed
541

542 543 544 545 546 547 548 549 550
#if AMR_MG_TIMER
    IpplTimings::stopTimer(bottomTimer_m);
#endif
        // phi = phi + e
        mglevel_m[level]->phi_p->update(1.0, *mglevel_m[level]->error_p, 1.0);
    }
}


551 552
void AmrMultiGrid::residual_no_fine_m(const lo_t& level,
                                      Teuchos::RCP<vector_t>& result,
553 554
                                      const Teuchos::RCP<vector_t>& rhs,
                                      const Teuchos::RCP<vector_t>& crs_rhs,
555
                                      const Teuchos::RCP<vector_t>& b)
556
{
frey_m's avatar
frey_m committed
557
    vector_t crse2fine(mglevel_m[level]->Anf_p->getDomainMap(), true);
kraus's avatar
kraus committed
558 559

    // get boundary for
560 561 562
    if ( mglevel_m[level]->Bcrse_p != Teuchos::null ) {
        mglevel_m[level]->Bcrse_p->apply(*crs_rhs, crse2fine);
    }
kraus's avatar
kraus committed
563

564 565 566 567 568
    // operation: crse2fine = 1.0 * crse2fine  + 1.0 * A^(l) * rhs
    mglevel_m[level]->Anf_p->apply(*rhs, crse2fine,
                                   Teuchos::NO_TRANS,
                                   scalar_t(1.0),
                                   scalar_t(1.0));
kraus's avatar
kraus committed
569

570
    result->update(1.0, *b, -1.0, crse2fine, 0.0);
571 572
}

573

frey_m's avatar
frey_m committed
574 575
#if AMR_MG_WRITE
void AmrMultiGrid::writeResidualNorm_m() {
576
    scalar_t err = 0.0;
kraus's avatar
kraus committed
577

578
    std::ofstream out;
frey_m's avatar
frey_m committed
579
    if ( comm_mp->getRank() == 0 )
580
        out.open("residual.dat", std::ios::app);
kraus's avatar
kraus committed
581

582
    for (int lev = 0; lev < nlevel_m; ++lev) {
583
        scalar_t tmp = evalNorm_m(mglevel_m[lev]->residual_p);
kraus's avatar
kraus committed
584

frey_m's avatar
frey_m committed
585
        if ( comm_mp->getRank() == 0 )
586 587
            out << std::setw(15) << std::right << tmp;
    }
kraus's avatar
kraus committed
588

frey_m's avatar
frey_m committed
589
    if ( comm_mp->getRank() == 0 )
590 591
        out.close();
}
frey_m's avatar
frey_m committed
592
#endif
593 594


595 596 597 598
typename AmrMultiGrid::scalar_t
AmrMultiGrid::evalNorm_m(const Teuchos::RCP<const vector_t>& x)
{
    scalar_t norm = 0.0;
kraus's avatar
kraus committed
599

600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616
    switch ( norm_m ) {
        case Norm::L1:
        {
            norm = x->norm1();
            break;
        }
        case Norm::L2:
        {
            norm = x->norm2();
            break;
        }
        case Norm::LINF:
        {
            norm = x->normInf();
            break;
        }
        default:
617 618
            throw OpalException("AmrMultiGrid::evalNorm_m()",
                                "This type of norm not suppported.");
619 620 621 622 623
    }
    return norm;
}


frey_m's avatar
frey_m committed
624 625 626 627 628
void AmrMultiGrid::initResidual_m(std::vector<scalar_t>& rhsNorms,
                                  std::vector<scalar_t>& resNorms)
{
    rhsNorms.clear();
    resNorms.clear();
kraus's avatar
kraus committed
629

630 631
#if AMR_MG_WRITE
    std::ofstream out;
kraus's avatar
kraus committed
632

frey_m's avatar
frey_m committed
633
    if ( comm_mp->getRank() == 0) {
634
        out.open("residual.dat", std::ios::out);
kraus's avatar
kraus committed
635

636 637 638 639 640
        for (int lev = 0; lev < nlevel_m; ++lev)
            out << std::setw(14) << std::right << "level" << lev;
        out << std::endl;
    }
#endif
kraus's avatar
kraus committed
641

642
    for (int lev = 0; lev < nlevel_m; ++lev) {
643 644
        this->residual_m(lev,
                         mglevel_m[lev]->residual_p,
645
                         mglevel_m[lev]->rho_p,
646
                         mglevel_m[lev]->phi_p);
kraus's avatar
kraus committed
647

frey_m's avatar
frey_m committed
648
        resNorms.push_back(evalNorm_m(mglevel_m[lev]->residual_p));
kraus's avatar
kraus committed
649

650
#if AMR_MG_WRITE
frey_m's avatar
frey_m committed
651
        if ( comm_mp->getRank() == 0 )
frey_m's avatar
frey_m committed
652
            out << std::setw(15) << std::right << resNorms.back();
653
#endif
kraus's avatar
kraus committed
654

frey_m's avatar
frey_m committed
655
        rhsNorms.push_back(evalNorm_m(mglevel_m[lev]->rho_p));
656
    }
kraus's avatar
kraus committed
657

658
#if AMR_MG_WRITE
frey_m's avatar
frey_m committed
659
    if ( comm_mp->getRank() == 0 )
660 661 662 663 664
        out.close();
#endif
}


665
void AmrMultiGrid::computeEfield_m(AmrVectorFieldContainer_t& efield) {
666 667 668
#if AMR_MG_TIMER
    IpplTimings::startTimer(efieldTimer_m);
#endif
669 670 671
    Teuchos::RCP<vector_t> efield_p = Teuchos::null;
    for (int lev = nlevel_m - 1; lev > -1; --lev) {
        int ilev = lbase_m + lev;
kraus's avatar
kraus committed
672

673
        efield_p = Teuchos::rcp( new vector_t(mglevel_m[lev]->map_p, false) );
kraus's avatar
kraus committed
674 675


676
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
677
            efield[ilev][d]->setVal(0.0, efield[ilev][d]->nGrow());
678
            efield_p->putScalar(0.0);
679
            mglevel_m[lev]->G_p[d]->apply(*mglevel_m[lev]->phi_p, *efield_p);
680
            this->trilinos2amrex_m(lev, 0, *efield[ilev][d], efield_p);
681 682
        }
    }
683 684 685
#if AMR_MG_TIMER
    IpplTimings::stopTimer(efieldTimer_m);
#endif
686 687 688
}


frey_m's avatar
frey_m committed
689 690
void AmrMultiGrid::setup_m(const amrex::Vector<AmrField_u>& rho,
                           const amrex::Vector<AmrField_u>& phi,
691
                           const bool& matrices)
692 693 694 695
{
#if AMR_MG_TIMER
    IpplTimings::startTimer(buildTimer_m);
#endif
frey_m's avatar
frey_m committed
696

697
    if ( lbase_m == lfine_m )
698
        this->buildSingleLevel_m(rho, phi, matrices);
699 700
    else
        this->buildMultiLevel_m(rho, phi, matrices);
kraus's avatar
kraus committed
701

702
    mglevel_m[lfine_m]->error_p->putScalar(0.0);
frey_m's avatar
frey_m committed
703

704
    if ( matrices ) {
frey_m's avatar
frey_m committed
705
        this->clearMasks_m();
706
        // set the bottom solve operator
707 708 709
        if (!solver_mp->hasOperator()) {
            solver_mp->setOperator(mglevel_m[lbase_m]->Anf_p, mglevel_m[0].get());
        }
710 711 712 713 714 715 716 717
    }

#if AMR_MG_TIMER
    IpplTimings::stopTimer(buildTimer_m);
#endif
}


frey_m's avatar
frey_m committed
718 719
void AmrMultiGrid::buildSingleLevel_m(const amrex::Vector<AmrField_u>& rho,
                                      const amrex::Vector<AmrField_u>& phi,
720 721 722
                                      const bool& matrices)
{
    this->open_m(lbase_m, matrices);
frey_m's avatar
frey_m committed
723

frey_m's avatar
frey_m committed
724
    const scalar_t* invdx = mglevel_m[lbase_m]->invCellSize();
kraus's avatar
kraus committed
725

726 727 728 729 730 731
    const scalar_t invdx2[] = {
        D_DECL( invdx[0] * invdx[0],
                invdx[1] * invdx[1],
                invdx[2] * invdx[2] )
    };

732
    if (matrices) {
733
        for (amrex::MFIter mfi(*mglevel_m[lbase_m]->mask, true);
734 735
             mfi.isValid(); ++mfi)
        {
736 737 738 739
            const box_t&       tbx = mfi.tilebox();
            const basefab_t&  mfab = (*mglevel_m[lbase_m]->mask)[mfi];
            const farraybox_t& rhofab = (*rho[lbase_m])[mfi];
            const farraybox_t& pfab = (*phi[lbase_m])[mfi];
kraus's avatar
kraus committed
740

741 742
            const int* lo = tbx.loVect();
            const int* hi = tbx.hiVect();
kraus's avatar
kraus committed
743

744 745 746 747 748 749
            for (int i = lo[0]; i <= hi[0]; ++i) {
                for (int j = lo[1]; j <= hi[1]; ++j) {
#if AMREX_SPACEDIM == 3
                    for (int k = lo[2]; k <= hi[2]; ++k) {
#endif
                        AmrIntVect_t iv(D_DECL(i, j, k));
frey_m's avatar
frey_m committed
750
                        go_t gidx = mglevel_m[lbase_m]->serialize(iv);
kraus's avatar
kraus committed
751

752 753 754 755
                        if (!solver_mp->hasOperator()) {
                            this->buildNoFinePoissonMatrix_m(lbase_m, gidx, iv, mfab, invdx2);
                            this->buildGradientMatrix_m(lbase_m, gidx, iv, mfab, invdx);
                        }
kraus's avatar
kraus committed
756

757 758
                        mglevel_m[lbase_m]->rho_p->replaceGlobalValue(gidx, rhofab(iv, 0));
                        mglevel_m[lbase_m]->phi_p->replaceGlobalValue(gidx, pfab(iv, 0));
kraus's avatar
kraus committed
759

760
#if AMREX_SPACEDIM == 3
761
                    }
762
#endif
763 764 765
                }
            }
        }
766 767 768 769 770
    } else {
        this->buildDensityVector_m(lbase_m, *rho[lbase_m]);
    }

    this->close_m(lbase_m, matrices);
frey_m's avatar
frey_m committed
771

772
    if (matrices) {
773
        mglevel_m[lbase_m]->Awf_p = Teuchos::null;
774 775
        mglevel_m[lbase_m]->UnCovered_p = Teuchos::null;
    }
776 777 778
}


frey_m's avatar
frey_m committed
779 780
void AmrMultiGrid::buildMultiLevel_m(const amrex::Vector<AmrField_u>& rho,
                                     const amrex::Vector<AmrField_u>& phi,
781 782 783
                                     const bool& matrices)
{
    // the base level has no smoother --> nlevel_m - 1
784 785 786 787 788
    if ( matrices ) {
        // although we do a resize afterwards, we do this to be safe
        for (int lev = nlevel_m-1; lev < (int)smoother_m.size(); ++lev) {
            smoother_m[lev].reset();
        }
789
        smoother_m.resize(nlevel_m-1);
790 791
    }
    
792
    for (int lev = 0; lev < nlevel_m; ++lev) {
793 794 795 796 797
        this->open_m(lev, matrices);

        int ilev = lbase_m + lev;

        // find all coarse cells that are covered by fine cells
frey_m's avatar
frey_m committed
798
//         AmrIntVect_t rr = mglevel_m[lev]->refinement();
799

frey_m's avatar
frey_m committed
800
        const scalar_t* invdx = mglevel_m[lev]->invCellSize();
801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817

        const scalar_t invdx2[] = {
            D_DECL( invdx[0] * invdx[0],
                    invdx[1] * invdx[1],
                    invdx[2] * invdx[2] )
        };

        if ( matrices ) {
            for (amrex::MFIter mfi(*mglevel_m[lev]->mask, true);
                 mfi.isValid(); ++mfi)
            {
                const box_t&       tbx = mfi.tilebox();
                const basefab_t&  mfab = (*mglevel_m[lev]->mask)[mfi];
                const basefab_t&  rfab = (*mglevel_m[lev]->refmask)[mfi];
                const basefab_t&  cfab = (*mglevel_m[lev]->crsemask)[mfi];
                const farraybox_t& rhofab = (*rho[ilev])[mfi];
                const farraybox_t& pfab = (*phi[ilev])[mfi];
kraus's avatar
kraus committed
818

819 820
                const int* lo = tbx.loVect();
                const int* hi = tbx.hiVect();
kraus's avatar
kraus committed
821

822 823 824 825 826 827 828 829 830
                for (int i = lo[0]; i <= hi[0]; ++i) {
                    int ii = i << 1;
                    for (int j = lo[1]; j <= hi[1]; ++j) {
                        int jj = j << 1;
#if AMREX_SPACEDIM == 3
                        for (int k = lo[2]; k <= hi[2]; ++k) {
                            int kk = k << 1;
#endif
                            AmrIntVect_t iv(D_DECL(i, j, k));