MGPoissonSolver.cpp 25 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9
//
// Class MGPoissonSolver
//   This class contains methods for solving Poisson's equation for the
//   space charge portion of the calculation.
//
//   A smoothed aggregation based AMG preconditioned iterative solver for space charge
//   \see FFTPoissonSolver
//   \warning This solver is in an EXPERIMENTAL STAGE. For reliable simulations use the FFTPoissonSolver
//
10 11 12
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
13 14
// All rights reserved
//
15 16 17 18 19
// Implemented as part of the master thesis
// "A Parallel Multigrid Solver for Beam Dynamics"
// and the paper
// "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations"
// (https://doi.org/10.1016/j.jcp.2010.02.022)
frey_m's avatar
frey_m committed
20
//
frey_m's avatar
frey_m committed
21
// In 2020, the code was ported to the second generation Trilinos packages,
frey_m's avatar
frey_m committed
22
// i.e., Epetra --> Tpetra, ML --> MueLu. See also issue #507.
frey_m's avatar
frey_m committed
23
//
frey_m's avatar
frey_m committed
24 25 26 27 28 29 30 31 32 33
// 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/>.
//
34

35
//#define DBG_STENCIL
36

kraus's avatar
kraus committed
37
#include "Solvers/MGPoissonSolver.h"
38 39 40 41 42 43

#include "Structure/BoundaryGeometry.h"
#include "ArbitraryDomain.h"
#include "EllipticDomain.h"
#include "BoxCornerDomain.h"

44
#include "Track/Track.h"
gsell's avatar
gsell committed
45 46 47 48 49
#include "Physics/Physics.h"
#include "Utilities/OpalException.h"
#include "Attributes/Attributes.h"
#include "ValueDefinitions/RealVariable.h"
#include "AbstractObjects/OpalData.h"
50
#include "Utilities/Options.h"
51

52 53
#include <Tpetra_Import.hpp>
#include <BelosTpetraAdapter.hpp>
54

frey_m's avatar
frey_m committed
55 56 57 58
#ifdef DBG_STENCIL
    #include "TpetraExt_MatrixMatrix.hpp"
#endif

59 60 61 62 63
#include "Teuchos_CommandLineProcessor.hpp"

#include "BelosLinearProblem.hpp"
#include "BelosRCGSolMgr.hpp"
#include "BelosBlockCGSolMgr.hpp"
frey_m's avatar
frey_m committed
64 65 66
#include "BelosBiCGStabSolMgr.hpp"
#include "BelosBlockGmresSolMgr.hpp"
#include "BelosGCRODRSolMgr.hpp"
67

68
#include <MueLu_CreateTpetraPreconditioner.hpp>
69

gsell's avatar
gsell committed
70 71
#include <algorithm>

72 73
using Teuchos::RCP;
using Teuchos::rcp;
gsell's avatar
gsell committed
74 75
using Physics::c;

76
MGPoissonSolver::MGPoissonSolver ( PartBunch *beam,
77 78 79 80 81 82 83
                                   Mesh_t *mesh,
                                   FieldLayout_t *fl,
                                   std::vector<BoundaryGeometry *> geometries,
                                   std::string itsolver,
                                   std::string interpl,
                                   double tol,
                                   int maxiters,
84
                                   std::string precmode)
frey_m's avatar
fix bug  
frey_m committed
85
    : isMatrixfilled_m(false)
86
    , useLeftPrec_m(true)
frey_m's avatar
fix bug  
frey_m committed
87
    , geometries_m(geometries)
88 89 90 91 92 93 94
    , tol_m(tol)
    , maxiters_m(maxiters)
    , comm_mp(new Comm_t(Ippl::getComm()))
    , itsBunch_m(beam)
    , mesh_m(mesh)
    , layout_m(fl)
{
95

gsell's avatar
gsell committed
96 97
    domain_m = layout_m->getDomain();

98
    for (int i = 0; i < 3; i++) {
gsell's avatar
gsell committed
99 100 101 102 103
        hr_m[i] = mesh_m->get_meshSpacing(i);
        orig_nr_m[i] = domain_m[i].length();
    }

    precmode_m = STD_PREC;
frey_m's avatar
frey_m committed
104
    if (precmode == "HIERARCHY") precmode_m = REUSE_HIERARCHY;
frey_m's avatar
frey_m committed
105
    else if (precmode == "REUSE") precmode_m = REUSE_PREC;
gsell's avatar
gsell committed
106

107
    repartFreq_m = 1000;
frey_m's avatar
frey_m committed
108
    if (Ippl::Info->getOutputLevel() > 3)
gsell's avatar
gsell committed
109 110 111 112
        verbose_m = true;
    else
        verbose_m = false;

113
    // Find CURRENT geometry
gsell's avatar
gsell committed
114
    currentGeometry = geometries_m[0];
frey_m's avatar
frey_m committed
115 116 117 118 119
    if (currentGeometry->getFilename() == "") {
        if (currentGeometry->getTopology() == "ELLIPTIC"){
            bp_m = std::unique_ptr<IrregularDomain>(
                new EllipticDomain(currentGeometry, orig_nr_m, hr_m, interpl));

kraus's avatar
kraus committed
120
        } else if (currentGeometry->getTopology() == "BOXCORNER") {
frey_m's avatar
frey_m committed
121 122 123 124 125 126 127 128
            bp_m = std::unique_ptr<IrregularDomain>(
                new BoxCornerDomain(currentGeometry->getA(),
                                    currentGeometry->getB(),
                                    currentGeometry->getC(),
                                    currentGeometry->getLength(),
                                    currentGeometry->getL1(),
                                    currentGeometry->getL2(),
                                    orig_nr_m, hr_m, interpl));
frey_m's avatar
frey_m committed
129
            bp_m->compute(itsBunch_m->get_hr());
gsell's avatar
gsell committed
130
        } else {
131 132
            throw OpalException("MGPoissonSolver::MGPoissonSolver",
                                "Geometry not known");
gsell's avatar
gsell committed
133
        }
134 135 136 137
    } else {
        NDIndex<3> localId = layout_m->getLocalNDIndex();
        if (localId[0].length() != domain_m[0].length() ||
            localId[1].length() != domain_m[1].length()) {
138
            throw OpalException("MGPoissonSolver::MGPoissonSolver",
139 140 141 142 143
                                "The class ArbitraryDomain only works with parallelization\n"
                                "in z-direction.\n"
                                "Please set PARFFTX=FALSE, PARFFTY=FALSE, PARFFTT=TRUE in \n"
                                "the definition of the field solver in the input file.\n");
        }
frey_m's avatar
frey_m committed
144 145
        bp_m = std::unique_ptr<IrregularDomain>(
            new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl));
146
    }
147

148
    map_p = Teuchos::null;
gsell's avatar
gsell committed
149 150 151
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
152
    prec_mp = Teuchos::null;
153

gsell's avatar
gsell committed
154 155
    numBlocks_m = Options::numBlocks;
    recycleBlocks_m = Options::recycleBlocks;
156
    nLHS_m = Options::nLHS;
frey_m's avatar
frey_m committed
157
    setupMueLuList();
frey_m's avatar
frey_m committed
158
    setupBelosList();
frey_m's avatar
frey_m committed
159 160 161
    problem_mp = rcp(new Belos::LinearProblem<TpetraScalar_t,
                                              TpetraMultiVector_t,
                                              TpetraOperator_t>);
162
    // setup Belos solver
frey_m's avatar
frey_m committed
163 164 165 166 167 168 169 170 171 172 173
    if (itsolver == "CG") {
        if (numBlocks_m == 0 || recycleBlocks_m == 0) {
            solver_mp = rcp(new Belos::BlockCGSolMgr<TpetraScalar_t,
                                                     TpetraMultiVector_t,
                                                     TpetraOperator_t>());
        } else {
            solver_mp = rcp(new Belos::RCGSolMgr<TpetraScalar_t,
                                                 TpetraMultiVector_t,
                                                 TpetraOperator_t>());
        }
    } else if (itsolver == "BICGSTAB") {
174
        useLeftPrec_m = false;
frey_m's avatar
frey_m committed
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
        solver_mp = rcp(new Belos::BiCGStabSolMgr<TpetraScalar_t,
                                                  TpetraMultiVector_t,
                                                  TpetraOperator_t>());
    } else if (itsolver == "GMRES") {
        if (numBlocks_m == 0 || recycleBlocks_m == 0) {
            solver_mp = rcp(new Belos::BlockGmresSolMgr<TpetraScalar_t,
                                                        TpetraMultiVector_t,
                                                        TpetraOperator_t>());
        } else {
            solver_mp = rcp(new Belos::GCRODRSolMgr<TpetraScalar_t,
                                                    TpetraMultiVector_t,
                                                    TpetraOperator_t>());
        }
    } else {
        throw OpalException("MGPoissonSolver", "No valid iterative solver selected!");
    }

192
    solver_mp->setParameters(rcp(&belosList, false));
frey_m's avatar
frey_m committed
193

gsell's avatar
gsell committed
194 195

    //all timers used here
196 197 198 199
    FunctionTimer1_m = IpplTimings::getTimer("BGF-IndexCoordMap");
    FunctionTimer2_m = IpplTimings::getTimer("computeMap");
    FunctionTimer3_m = IpplTimings::getTimer("IPPL to RHS");
    FunctionTimer4_m = IpplTimings::getTimer("ComputeStencil");
frey_m's avatar
frey_m committed
200
    FunctionTimer5_m = IpplTimings::getTimer("MueLu");
201 202 203
    FunctionTimer6_m = IpplTimings::getTimer("Setup");
    FunctionTimer7_m = IpplTimings::getTimer("CG");
    FunctionTimer8_m = IpplTimings::getTimer("LHS to IPPL");
gsell's avatar
gsell committed
204 205
}

206
void MGPoissonSolver::deletePtr() {
207
    map_p  = Teuchos::null;
208 209 210
    A      = Teuchos::null;
    LHS    = Teuchos::null;
    RHS    = Teuchos::null;
211
    prec_mp = Teuchos::null;
frey_m's avatar
frey_m committed
212
    isMatrixfilled_m = false;
213 214
}

215 216
MGPoissonSolver::~MGPoissonSolver() {
    deletePtr ();
217 218
    solver_mp = Teuchos::null;
    problem_mp = Teuchos::null;
gsell's avatar
gsell committed
219 220
}

221
void MGPoissonSolver::computePotential(Field_t& /*rho*/, Vector_t /*hr*/, double /*zshift*/) {
gsell's avatar
gsell committed
222 223 224
    throw OpalException("MGPoissonSolver", "not implemented yet");
}

225
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
226
    if (itsBunch_m->getLocalTrackStep()%repartFreq_m == 0){
227
        deletePtr();
frey_m's avatar
frey_m committed
228
        IPPLToMap3D(localId);
229 230

        extrapolateLHS();
231 232
    }
}
gsell's avatar
gsell committed
233

234 235 236 237 238
void MGPoissonSolver::extrapolateLHS() {
// Aitken-Neville
// Pi0 (x) := yi , i = 0 : n
// Pik (x) := (x − xi ) Pi+1,k−1(x) − (x − xi+k ) Pi,k−1(x) /(xi+k − xi )
// k = 1, . . . , n, i = 0, . . . , n − k.
239
    //we also have to redistribute LHS
240

241
    if (Teuchos::is_null(LHS)){
242 243
        LHS = rcp(new TpetraVector_t(map_p));
        LHS->putScalar(1.0);
244
    } else {
245 246 247
        RCP<TpetraVector_t> tmplhs = rcp(new TpetraVector_t(map_p));
        Tpetra::Import<> importer(map_p, LHS->getMap());
        tmplhs->doImport(*LHS, importer, Tpetra::CombineMode::ADD);
248 249
        LHS = tmplhs;
    }
gsell's avatar
gsell committed
250

251
    //...and all previously saved LHS
252
    std::deque< TpetraVector_t >::iterator it = OldLHS.begin();
frey_m's avatar
frey_m committed
253
    if (OldLHS.size() > 0) {
254
        int n = OldLHS.size();
frey_m's avatar
frey_m committed
255
        for (int i = 0; i < n; ++i) {
256 257 258
            TpetraVector_t tmplhs = TpetraVector_t(map_p);
            Tpetra::Import<> importer(map_p, it->getMap());
            tmplhs.doImport(*it, importer, Tpetra::CombineMode::ADD);
259 260 261 262
            *it = tmplhs;
            ++it;
        }
    }
gsell's avatar
gsell committed
263

264 265
    // extrapolate last OldLHS.size LHS to get a new start vector
    it = OldLHS.begin();
frey_m's avatar
frey_m committed
266
    if (nLHS_m == 0 || OldLHS.size()==0)
267
        LHS->putScalar(1.0);
frey_m's avatar
frey_m committed
268
    else if (OldLHS.size() == 1)
269
        *LHS = *it;
frey_m's avatar
frey_m committed
270
    else if (OldLHS.size() == 2)
271
        LHS->update(2.0, *it++, -1.0, *it, 0.0);
frey_m's avatar
frey_m committed
272
    else if (OldLHS.size() > 2) {
273
        int n = OldLHS.size();
274
        P_mp = rcp(new TpetraMultiVector_t(map_p, nLHS_m, false));
frey_m's avatar
frey_m committed
275
        for (int i = 0; i < n; ++i) {
276
           *(P_mp->getVectorNonConst(i)) = *it++;
277
        }
frey_m's avatar
frey_m committed
278 279
        for (int k = 1; k < n; ++k) {
           for (int i = 0; i < n - k; ++i) {
280
              P_mp->getVectorNonConst(i)->update(-(i + 1) / (float)k, *(P_mp->getVector(i + 1)), (i + k + 1) / (float)k);
281
           }
282
        }
283
        *LHS = *(P_mp->getVector(0));
284
     } else
frey_m's avatar
frey_m committed
285 286
        throw OpalException("MGPoissonSolver",
                            "Invalid number of old LHS: " + std::to_string(OldLHS.size()));
gsell's avatar
gsell committed
287 288
}

289

gsell's avatar
gsell committed
290 291
// given a charge-density field rho and a set of mesh spacings hr,
// compute the electric field and put in eg by solving the Poisson's equation
frey_m's avatar
frey_m committed
292
// XXX: use matrix stencil in computation directly (no Tpetra, define operators
gsell's avatar
gsell committed
293 294
// on IPPL GRID)
void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
295

kraus's avatar
kraus committed
296
    Inform msg("OPAL ", INFORM_ALL_NODES);
gsell's avatar
gsell committed
297 298 299
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
300

frey_m's avatar
frey_m committed
301 302 303
    bp_m->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
    bp_m->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
    bp_m->setNr(nr_m);
304

305 306
    NDIndex<3> localId = layout_m->getLocalNDIndex();

307
    IpplTimings::startTimer(FunctionTimer1_m);
308
    // Compute boundary intersections (only do on the first step)
frey_m's avatar
frey_m committed
309
    if (!itsBunch_m->getLocalTrackStep())
frey_m's avatar
frey_m committed
310
        bp_m->compute(hr, localId);
311
    IpplTimings::stopTimer(FunctionTimer1_m);
312 313

    // Define the Map
kraus's avatar
kraus committed
314
    INFOMSG(level3 << "* Computing Map..." << endl);
315
    IpplTimings::startTimer(FunctionTimer2_m);
316
    computeMap(localId);
317
    IpplTimings::stopTimer(FunctionTimer2_m);
kraus's avatar
kraus committed
318
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
319

frey_m's avatar
frey_m committed
320
    // Allocate the RHS with the new Tpetra Map
321
    if (Teuchos::is_null(RHS))
322 323
        RHS = rcp(new TpetraVector_t(map_p));
    RHS->putScalar(0.0);
324

frey_m's avatar
frey_m committed
325
    // get charge densities from IPPL field and store in Tpetra vector (RHS)
kraus's avatar
kraus committed
326 327 328 329 330
    if (verbose_m) {
        Ippl::Comm->barrier();
        msg << "* Node:" << Ippl::myNode() << ", Filling RHS..." << endl;
        Ippl::Comm->barrier();
    }
331 332
    IpplTimings::startTimer(FunctionTimer3_m);
    int id = 0;
333
    float scaleFactor = itsBunch_m->getdT();
334

335

kraus's avatar
kraus committed
336
    if (verbose_m) {
frey_m's avatar
frey_m committed
337 338 339
        msg << "* Node:" << Ippl::myNode() << ", Rho for final element: "
            << rho[localId[0].last()][localId[1].last()][localId[2].last()].get()
            << endl;
340

kraus's avatar
kraus committed
341
        Ippl::Comm->barrier();
frey_m's avatar
frey_m committed
342 343 344 345 346
        msg << "* Node:" << Ippl::myNode() << ", Local nx*ny*nz = "
            <<  localId[2].last() *  localId[0].last() *  localId[1].last()
            << endl;
        msg << "* Node:" << Ippl::myNode()
            << ", Number of reserved local elements in RHS: "
347
            << RHS->getLocalLength() << endl;
frey_m's avatar
frey_m committed
348 349
        msg << "* Node:" << Ippl::myNode()
            << ", Number of reserved global elements in RHS: "
350
            << RHS->getGlobalLength() << endl;
kraus's avatar
kraus committed
351 352
        Ippl::Comm->barrier();
    }
353
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
354
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
kraus's avatar
kraus committed
355
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
356
                if (bp_m->isInside(idx, idy, idz))
357
                        RHS->getDataNonConst()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
gsell's avatar
gsell committed
358 359 360
            }
        }
    }
361

362
    IpplTimings::stopTimer(FunctionTimer3_m);
kraus's avatar
kraus committed
363 364
    if (verbose_m) {
        Ippl::Comm->barrier();
frey_m's avatar
frey_m committed
365 366
        msg << "* Node:" << Ippl::myNode()
            << ", Number of Local Inside Points " << id << endl;
kraus's avatar
kraus committed
367 368 369
        msg << "* Node:" << Ippl::myNode() << ", Done." << endl;
        Ippl::Comm->barrier();
    }
gsell's avatar
gsell committed
370
    // build discretization matrix
kraus's avatar
kraus committed
371
    INFOMSG(level3 << "* Building Discretization Matrix..." << endl);
372
    IpplTimings::startTimer(FunctionTimer4_m);
frey_m's avatar
frey_m committed
373
    if (Teuchos::is_null(A))
374
        A = rcp(new TpetraCrsMatrix_t(map_p,  7, Tpetra::StaticProfile));
375
    ComputeStencil(hr, RHS);
376
    IpplTimings::stopTimer(FunctionTimer4_m);
kraus's avatar
kraus committed
377
    INFOMSG(level3 << "* Done." << endl);
378

gsell's avatar
gsell committed
379
#ifdef DBG_STENCIL
frey_m's avatar
frey_m committed
380 381
    Tpetra::MatrixMarket::Writer<TpetraCrsMatrix_t>::writeSparseFile(
        "DiscrStencil.dat", A);
gsell's avatar
gsell committed
382
#endif
383

kraus's avatar
kraus committed
384
    INFOMSG(level3 << "* Computing Preconditioner..." << endl);
385
    IpplTimings::startTimer(FunctionTimer5_m);
frey_m's avatar
frey_m committed
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402
    if (Teuchos::is_null(prec_mp)) {
        Teuchos::RCP<TpetraOperator_t> At = Teuchos::rcp_dynamic_cast<TpetraOperator_t>(A);
        prec_mp = MueLu::CreateTpetraPreconditioner(At, MueLuList_m);
    }

    switch (precmode_m) {
        case REUSE_PREC:
        case REUSE_HIERARCHY: {
            MueLu::ReuseTpetraPreconditioner(A, *prec_mp);
            break;
        }
        case STD_PREC:
        default: {
            Teuchos::RCP<TpetraOperator_t> At = Teuchos::rcp_dynamic_cast<TpetraOperator_t>(A);
            prec_mp = MueLu::CreateTpetraPreconditioner(At, MueLuList_m);
            break;
        }
403
    }
404
    IpplTimings::stopTimer(FunctionTimer5_m);
kraus's avatar
kraus committed
405
    INFOMSG(level3 << "* Done." << endl);
406

gsell's avatar
gsell committed
407 408
    // setup preconditioned iterative solver
    // use old LHS solution as initial guess
kraus's avatar
kraus committed
409
    INFOMSG(level3 << "* Final Setup of Solver..." << endl);
410
    IpplTimings::startTimer(FunctionTimer6_m);
411 412 413
    problem_mp->setOperator(A);
    problem_mp->setLHS(LHS);
    problem_mp->setRHS(RHS);
414 415 416 417 418 419

    if (useLeftPrec_m)
        problem_mp->setLeftPrec(prec_mp);
    else
        problem_mp->setRightPrec(prec_mp);

420 421 422
    solver_mp->setProblem( problem_mp);
    if (!problem_mp->isProblemSet()) {
        if (problem_mp->setProblem() == false) {
423
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
424
        }
gsell's avatar
gsell committed
425
    }
426
    IpplTimings::stopTimer(FunctionTimer6_m);
kraus's avatar
kraus committed
427
    INFOMSG(level3 << "* Done." << endl);
428

gsell's avatar
gsell committed
429 430
    double time = MPI_Wtime();

kraus's avatar
kraus committed
431
    INFOMSG(level3 << "* Solving for Space Charge..." << endl);
432
    IpplTimings::startTimer(FunctionTimer7_m);
433
    solver_mp->solve();
frey_m's avatar
frey_m committed
434

435
    IpplTimings::stopTimer(FunctionTimer7_m);
kraus's avatar
kraus committed
436
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
437

kraus's avatar
kraus committed
438 439 440
    std::ofstream timings;
    if (true || verbose_m) {
        time = MPI_Wtime() - time;
441
        double minTime = 0, maxTime = 0, avgTime = 0;
442 443 444
        reduce(time, minTime, 1, std::less<double>());
        reduce(time, maxTime, 1, std::greater<double>());
        reduce(time, avgTime, 1, std::plus<double>());
445 446
        avgTime /= comm_mp->getSize();
        if (comm_mp->getRank() == 0) {
kraus's avatar
kraus committed
447
            char filename[50];
frey_m's avatar
frey_m committed
448 449
            sprintf(filename, "timing_MX%d_MY%d_MZ%d_nProc%d_recB%d_numB%d_nLHS%d",
                    orig_nr_m[0], orig_nr_m[1], orig_nr_m[2],
450
                    comm_mp->getSize(), recycleBlocks_m, numBlocks_m, nLHS_m);
frey_m's avatar
frey_m committed
451

kraus's avatar
kraus committed
452
            timings.open(filename, std::ios::app);
453
            timings << solver_mp->getNumIters() << "\t"
kraus's avatar
kraus committed
454 455 456 457 458 459 460 461 462 463 464 465 466 467
                    //<< time <<" "<<
                    << minTime << "\t"
                    << maxTime << "\t"
                    << avgTime << "\t"
                    << numBlocks_m << "\t"
                    << recycleBlocks_m << "\t"
                    << nLHS_m << "\t"
                    //<< OldLHS.size() <<"\t"<<
                    << std::endl;

            timings.close();
        }

    }
gsell's avatar
gsell committed
468
    // Store new LHS in OldLHS
frey_m's avatar
frey_m committed
469 470
    if (nLHS_m > 1) OldLHS.push_front(*(LHS.get()));
    if (OldLHS.size() > nLHS_m) OldLHS.pop_back();
gsell's avatar
gsell committed
471 472

    //now transfer solution back to IPPL grid
473 474
    IpplTimings::startTimer(FunctionTimer8_m);
    id = 0;
gsell's avatar
gsell committed
475
    rho = 0.0;
476
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
477
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
kraus's avatar
kraus committed
478
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
479
                  NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
frey_m's avatar
frey_m committed
480
                  if (bp_m->isInside(idx, idy, idz))
481
                     rho.localElement(l) = LHS->getData()[id++] * scaleFactor;
gsell's avatar
gsell committed
482 483 484
            }
        }
    }
485
    IpplTimings::stopTimer(FunctionTimer8_m);
486

frey_m's avatar
frey_m committed
487
    if (itsBunch_m->getLocalTrackStep()+1 == (long long)Track::block->localTimeSteps.front()) {
488 489 490
        A = Teuchos::null;
        LHS = Teuchos::null;
        RHS = Teuchos::null;
491
        prec_mp = Teuchos::null;
492
    }
gsell's avatar
gsell committed
493 494 495
}


496
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
497

gsell's avatar
gsell committed
498
    int NumMyElements = 0;
frey_m's avatar
frey_m committed
499
    std::vector<TpetraGlobalOrdinal_t> MyGlobalElements;
gsell's avatar
gsell committed
500

501
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
502
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
503
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
504 505
                if (bp_m->isInside(idx, idy, idz)) {
                    MyGlobalElements.push_back(bp_m->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
506 507 508
                    NumMyElements++;
                }
            }
kraus's avatar
kraus committed
509
        }
510
    }
511 512 513
    int indexbase = 0;
    map_p = Teuchos::rcp(new TpetraMap_t(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
                                         &MyGlobalElements[0], NumMyElements, indexbase, comm_mp));
gsell's avatar
gsell committed
514 515
}

516
void MGPoissonSolver::ComputeStencil(Vector_t /*hr*/, Teuchos::RCP<TpetraVector_t> RHS) {
gsell's avatar
gsell committed
517

frey_m's avatar
frey_m committed
518
    A->resumeFill();
519
    A->setAllToScalar(0.0);
520

521 522
    int NumMyElements = map_p->getNodeNumElements();
    auto MyGlobalElements = map_p->getMyGlobalIndices();
gsell's avatar
gsell committed
523

frey_m's avatar
frey_m committed
524 525
    std::vector<TpetraScalar_t> Values(6);
    std::vector<TpetraGlobalOrdinal_t> Indices(6);
gsell's avatar
gsell committed
526

frey_m's avatar
frey_m committed
527
    for (int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
528 529 530

        int NumEntries = 0;

531
        double WV, EV, SV, NV, FV, BV, CV, scaleFactor=1.0;
532
        int W, E, S, N, F, B;
gsell's avatar
gsell committed
533

frey_m's avatar
frey_m committed
534
        bp_m->getBoundaryStencil(MyGlobalElements[i], WV, EV, SV, NV, FV, BV, CV, scaleFactor);
535
        RHS->scale(scaleFactor);
gsell's avatar
gsell committed
536

frey_m's avatar
frey_m committed
537
        bp_m->getNeighbours(MyGlobalElements[i], W, E, S, N, F, B);
frey_m's avatar
frey_m committed
538
        if (E != -1) {
gsell's avatar
gsell committed
539 540 541
            Indices[NumEntries] = E;
            Values[NumEntries++] = EV;
        }
frey_m's avatar
frey_m committed
542
        if (W != -1) {
gsell's avatar
gsell committed
543 544 545
            Indices[NumEntries] = W;
            Values[NumEntries++] = WV;
        }
frey_m's avatar
frey_m committed
546
        if (S != -1) {
gsell's avatar
gsell committed
547 548 549
            Indices[NumEntries] = S;
            Values[NumEntries++] = SV;
        }
frey_m's avatar
frey_m committed
550
        if (N != -1) {
gsell's avatar
gsell committed
551 552 553
            Indices[NumEntries] = N;
            Values[NumEntries++] = NV;
        }
frey_m's avatar
frey_m committed
554
        if (F != -1) {
gsell's avatar
gsell committed
555 556 557
            Indices[NumEntries] = F;
            Values[NumEntries++] = FV;
        }
frey_m's avatar
frey_m committed
558
        if (B != -1) {
gsell's avatar
gsell committed
559 560 561 562
            Indices[NumEntries] = B;
            Values[NumEntries++] = BV;
        }

563
        // if matrix has already been filled (fillComplete()) we can only
gsell's avatar
gsell committed
564
        // replace entries
565

frey_m's avatar
fix bug  
frey_m committed
566
        if (isMatrixfilled_m) {
567
            // off-diagonal entries
568
            A->replaceGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
569
            // diagonal entry
570
            A->replaceGlobalValues(MyGlobalElements[i], 1, &CV, &MyGlobalElements[i]);
571
        } else {
gsell's avatar
gsell committed
572
            // off-diagonal entries
573
            A->insertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
gsell's avatar
gsell committed
574
            // diagonal entry
575
            A->insertGlobalValues(MyGlobalElements[i], 1, &CV, &MyGlobalElements[i]);
576
        }
gsell's avatar
gsell committed
577 578
    }

579 580 581
    RCP<ParameterList_t> params = Teuchos::parameterList();
    params->set ("Optimize Storage", true);
    A->fillComplete(params);
frey_m's avatar
fix bug  
frey_m committed
582
    isMatrixfilled_m = true;
gsell's avatar
gsell committed
583 584 585 586 587
}

void MGPoissonSolver::printLoadBalanceStats() {

    //compute some load balance statistics
588 589
    size_t myNumPart = map_p->getNodeNumElements();
    size_t NumPart = map_p->getGlobalNumElements() * 1.0 / comm_mp->getSize();
gsell's avatar
gsell committed
590
    double imbalance = 1.0;
frey_m's avatar
frey_m committed
591
    if (myNumPart >= NumPart)
gsell's avatar
gsell committed
592 593 594 595 596
        imbalance += (myNumPart - NumPart) / NumPart;
    else
        imbalance += (NumPart - myNumPart) / NumPart;

    double max = 0.0, min = 0.0, avg = 0.0;
frey_m's avatar
frey_m committed
597 598 599 600 601 602
    size_t minn = 0, maxn = 0;
    reduce(imbalance, min, 1, std::less<double>());
    reduce(imbalance, max, 1, std::greater<double>());
    reduce(imbalance, avg, 1, std::plus<double>());
    reduce(myNumPart, minn, 1, std::less<size_t>());
    reduce(myNumPart, maxn, 1, std::greater<size_t>());
gsell's avatar
gsell committed
603

604
    avg /= comm_mp->getSize();
kraus's avatar
kraus committed
605 606
    *gmsg << "LBAL min = " << min << ", max = " << max << ", avg = " << avg << endl;
    *gmsg << "min nr gridpoints = " << minn << ", max nr gridpoints = " << maxn << endl;
gsell's avatar
gsell committed
607 608
}

frey_m's avatar
frey_m committed
609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 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
void MGPoissonSolver::setupBelosList() {
    belosList.set("Maximum Iterations", maxiters_m);
    belosList.set("Convergence Tolerance", tol_m);

    if (numBlocks_m != 0 && recycleBlocks_m != 0){//only set if solver==RCGSolMgr
        belosList.set("Num Blocks", numBlocks_m);               // Maximum number of blocks in Krylov space
        belosList.set("Num Recycled Blocks", recycleBlocks_m); // Number of vectors in recycle space
    }
    if (verbose_m) {
        belosList.set("Verbosity", Belos::Errors + Belos::Warnings +
                                   Belos::TimingDetails + Belos::FinalSummary +
                                   Belos::StatusTestDetails);
        belosList.set("Output Frequency", 1);
    } else
        belosList.set("Verbosity", Belos::Errors + Belos::Warnings);
}


void MGPoissonSolver::setupMueLuList() {
    MueLuList_m.set("problem: type", "Poisson-3D");
    MueLuList_m.set("verbosity", "none");
    MueLuList_m.set("number of equations", 1);
    MueLuList_m.set("max levels", 8);
    MueLuList_m.set("cycle type", "V");

    // heuristic for max coarse size depending on number of processors
    int coarsest_size = std::max(comm_mp->getSize() * 10, 1024);
    MueLuList_m.set("coarse: max size", coarsest_size);

    MueLuList_m.set("multigrid algorithm", "sa");
    MueLuList_m.set("sa: damping factor", 1.33);
    MueLuList_m.set("sa: use filtered matrix", true);
    MueLuList_m.set("filtered matrix: reuse eigenvalue", false);

    MueLuList_m.set("repartition: enable", false);
    MueLuList_m.set("repartition: rebalance P and R", false);
    MueLuList_m.set("repartition: partitioner", "zoltan2");
    MueLuList_m.set("repartition: min rows per proc", 800);
    MueLuList_m.set("repartition: start level", 2);

    MueLuList_m.set("smoother: type", "CHEBYSHEV");
    MueLuList_m.set("smoother: pre or post", "both");
    Teuchos::ParameterList smparms;
    smparms.set("chebyshev: degree", 3);
    smparms.set("chebyshev: assume matrix does not change", false);
    smparms.set("chebyshev: zero starting solution", true);
    smparms.set("relaxation: sweeps", 3);
    MueLuList_m.set("smoother: params", smparms);

    MueLuList_m.set("smoother: type", "CHEBYSHEV");
    MueLuList_m.set("smoother: pre or post", "both");

    MueLuList_m.set("coarse: type", "KLU2");

    MueLuList_m.set("aggregation: type", "uncoupled");
    MueLuList_m.set("aggregation: min agg size", 3);
    MueLuList_m.set("aggregation: max agg size", 27);

    MueLuList_m.set("transpose: use implicit", false);

    switch (precmode_m) {
        case REUSE_PREC:
            MueLuList_m.set("reuse: type", "full");
            break;
        case REUSE_HIERARCHY:
            MueLuList_m.set("sa: use filtered matrix", false);
            MueLuList_m.set("reuse: type", "tP");
            break;
        case STD_PREC:
        default:
            MueLuList_m.set("reuse: type", "none");
            break;
    }
}

gsell's avatar
gsell committed
684 685 686 687
Inform &MGPoissonSolver::print(Inform &os) const {
    os << "* *************** M G P o i s s o n S o l v e r ************************************ " << endl;
    os << "* h " << hr_m << '\n';
    os << "* ********************************************************************************** " << endl;
688
    return os;
gsell's avatar
gsell committed
689 690
}

691 692 693 694 695 696 697
// vi: set et ts=4 sw=4 sts=4:
// Local Variables:
// mode:c
// c-basic-offset: 4
// indent-tabs-mode: nil
// require-final-newline: nil
// End: