MGPoissonSolver.cpp 25.4 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

#include "Structure/BoundaryGeometry.h"
#include "ArbitraryDomain.h"
#include "EllipticDomain.h"
#include "BoxCornerDomain.h"
43
#include "RectangularDomain.h"
44

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

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

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

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

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

69
#include <MueLu_CreateTpetraPreconditioner.hpp>
70

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

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

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

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

99
    for (int i = 0; i < 3; i++) {
gsell's avatar
gsell committed
100 101 102 103 104
        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
105
    if (precmode == "HIERARCHY") precmode_m = REUSE_HIERARCHY;
frey_m's avatar
frey_m committed
106
    else if (precmode == "REUSE") precmode_m = REUSE_PREC;
gsell's avatar
gsell committed
107

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

114
    // Find CURRENT geometry
gsell's avatar
gsell committed
115
    currentGeometry = geometries_m[0];
frey_m's avatar
frey_m committed
116 117 118 119 120
    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
121
        } else if (currentGeometry->getTopology() == "BOXCORNER") {
frey_m's avatar
frey_m committed
122 123 124 125 126 127 128
            bp_m = std::unique_ptr<IrregularDomain>(
                new BoxCornerDomain(currentGeometry->getA(),
                                    currentGeometry->getB(),
                                    currentGeometry->getC(),
                                    currentGeometry->getL1(),
                                    currentGeometry->getL2(),
                                    orig_nr_m, hr_m, interpl));
129
            bp_m->compute(itsBunch_m->get_hr(), layout_m->getLocalNDIndex());
frey_m's avatar
frey_m committed
130
        } else if (currentGeometry->getTopology() == "RECTANGULAR") {
131 132 133 134
            bp_m = std::unique_ptr<IrregularDomain>(
                new RectangularDomain(currentGeometry->getA(),
                                      currentGeometry->getB(),
                                      orig_nr_m, hr_m));
gsell's avatar
gsell committed
135
        } else {
136 137
            throw OpalException("MGPoissonSolver::MGPoissonSolver",
                                "Geometry not known");
gsell's avatar
gsell committed
138
        }
139 140 141 142
    } else {
        NDIndex<3> localId = layout_m->getLocalNDIndex();
        if (localId[0].length() != domain_m[0].length() ||
            localId[1].length() != domain_m[1].length()) {
143
            throw OpalException("MGPoissonSolver::MGPoissonSolver",
144 145 146 147 148
                                "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
149 150
        Vector_t hr = (currentGeometry->getmaxcoords() -
                       currentGeometry->getmincoords()) / orig_nr_m;
frey_m's avatar
frey_m committed
151
        bp_m = std::unique_ptr<IrregularDomain>(
frey_m's avatar
frey_m committed
152
            new ArbitraryDomain(currentGeometry, orig_nr_m, hr, interpl));
153
    }
154

155
    map_p = Teuchos::null;
gsell's avatar
gsell committed
156 157 158
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
159
    prec_mp = Teuchos::null;
160

gsell's avatar
gsell committed
161 162
    numBlocks_m = Options::numBlocks;
    recycleBlocks_m = Options::recycleBlocks;
163
    nLHS_m = Options::nLHS;
frey_m's avatar
frey_m committed
164
    setupMueLuList();
frey_m's avatar
frey_m committed
165
    setupBelosList();
frey_m's avatar
frey_m committed
166 167 168
    problem_mp = rcp(new Belos::LinearProblem<TpetraScalar_t,
                                              TpetraMultiVector_t,
                                              TpetraOperator_t>);
169
    // setup Belos solver
frey_m's avatar
frey_m committed
170 171 172 173 174 175 176 177 178 179 180
    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") {
181
        useLeftPrec_m = false;
frey_m's avatar
frey_m committed
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
        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!");
    }

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

gsell's avatar
gsell committed
201 202

    //all timers used here
203 204 205 206
    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
207
    FunctionTimer5_m = IpplTimings::getTimer("MueLu");
208 209 210
    FunctionTimer6_m = IpplTimings::getTimer("Setup");
    FunctionTimer7_m = IpplTimings::getTimer("CG");
    FunctionTimer8_m = IpplTimings::getTimer("LHS to IPPL");
gsell's avatar
gsell committed
211 212
}

213
void MGPoissonSolver::deletePtr() {
214
    map_p  = Teuchos::null;
215 216 217
    A      = Teuchos::null;
    LHS    = Teuchos::null;
    RHS    = Teuchos::null;
218
    prec_mp = Teuchos::null;
frey_m's avatar
frey_m committed
219
    isMatrixfilled_m = false;
220 221
}

222 223
MGPoissonSolver::~MGPoissonSolver() {
    deletePtr ();
224 225
    solver_mp = Teuchos::null;
    problem_mp = Teuchos::null;
gsell's avatar
gsell committed
226 227
}

228
void MGPoissonSolver::computePotential(Field_t& /*rho*/, Vector_t /*hr*/, double /*zshift*/) {
gsell's avatar
gsell committed
229 230 231
    throw OpalException("MGPoissonSolver", "not implemented yet");
}

232
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
233
    if (itsBunch_m->getLocalTrackStep()%repartFreq_m == 0){
234
        deletePtr();
frey_m's avatar
frey_m committed
235
        IPPLToMap3D(localId);
236 237

        extrapolateLHS();
238 239
    }
}
gsell's avatar
gsell committed
240

241 242 243 244 245
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.
246
    //we also have to redistribute LHS
247

248
    if (Teuchos::is_null(LHS)){
249 250
        LHS = rcp(new TpetraVector_t(map_p));
        LHS->putScalar(1.0);
251
    } else {
252 253 254
        RCP<TpetraVector_t> tmplhs = rcp(new TpetraVector_t(map_p));
        Tpetra::Import<> importer(map_p, LHS->getMap());
        tmplhs->doImport(*LHS, importer, Tpetra::CombineMode::ADD);
255 256
        LHS = tmplhs;
    }
gsell's avatar
gsell committed
257

258
    //...and all previously saved LHS
259
    std::deque< TpetraVector_t >::iterator it = OldLHS.begin();
frey_m's avatar
frey_m committed
260
    if (OldLHS.size() > 0) {
261
        int n = OldLHS.size();
frey_m's avatar
frey_m committed
262
        for (int i = 0; i < n; ++i) {
263 264 265
            TpetraVector_t tmplhs = TpetraVector_t(map_p);
            Tpetra::Import<> importer(map_p, it->getMap());
            tmplhs.doImport(*it, importer, Tpetra::CombineMode::ADD);
266 267 268 269
            *it = tmplhs;
            ++it;
        }
    }
gsell's avatar
gsell committed
270

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

296

gsell's avatar
gsell committed
297 298
// 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
299
// XXX: use matrix stencil in computation directly (no Tpetra, define operators
gsell's avatar
gsell committed
300 301
// on IPPL GRID)
void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
302

kraus's avatar
kraus committed
303
    Inform msg("OPAL ", INFORM_ALL_NODES);
gsell's avatar
gsell committed
304 305 306
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
307

frey_m's avatar
frey_m committed
308
    bp_m->setNr(nr_m);
309

310 311
    NDIndex<3> localId = layout_m->getLocalNDIndex();

312
    IpplTimings::startTimer(FunctionTimer1_m);
313
    // Compute boundary intersections (only do on the first step)
frey_m's avatar
frey_m committed
314
    if (!itsBunch_m->getLocalTrackStep())
frey_m's avatar
frey_m committed
315
        bp_m->compute(hr, localId);
316
    IpplTimings::stopTimer(FunctionTimer1_m);
317 318

    // Define the Map
kraus's avatar
kraus committed
319
    INFOMSG(level3 << "* Computing Map..." << endl);
320
    IpplTimings::startTimer(FunctionTimer2_m);
321
    computeMap(localId);
322
    IpplTimings::stopTimer(FunctionTimer2_m);
kraus's avatar
kraus committed
323
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
324

frey_m's avatar
frey_m committed
325
    // Allocate the RHS with the new Tpetra Map
326
    if (Teuchos::is_null(RHS))
327 328
        RHS = rcp(new TpetraVector_t(map_p));
    RHS->putScalar(0.0);
329

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

340

kraus's avatar
kraus committed
341
    if (verbose_m) {
frey_m's avatar
frey_m committed
342 343 344
        msg << "* Node:" << Ippl::myNode() << ", Rho for final element: "
            << rho[localId[0].last()][localId[1].last()][localId[2].last()].get()
            << endl;
345

kraus's avatar
kraus committed
346
        Ippl::Comm->barrier();
frey_m's avatar
frey_m committed
347 348 349 350 351
        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: "
352
            << RHS->getLocalLength() << endl;
frey_m's avatar
frey_m committed
353 354
        msg << "* Node:" << Ippl::myNode()
            << ", Number of reserved global elements in RHS: "
355
            << RHS->getGlobalLength() << endl;
kraus's avatar
kraus committed
356 357
        Ippl::Comm->barrier();
    }
358
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
359
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
kraus's avatar
kraus committed
360
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
361
                NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
frey_m's avatar
frey_m committed
362
                if (bp_m->isInside(idx, idy, idz))
frey_m's avatar
frey_m committed
363 364
                        RHS->replaceGlobalValue(bp_m->getIdx(idx, idy, idz),
                                                4.0 * M_PI * rho.localElement(l) / scaleFactor);
gsell's avatar
gsell committed
365 366 367
            }
        }
    }
368

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

gsell's avatar
gsell committed
386
#ifdef DBG_STENCIL
frey_m's avatar
frey_m committed
387 388
    Tpetra::MatrixMarket::Writer<TpetraCrsMatrix_t>::writeSparseFile(
        "DiscrStencil.dat", A);
gsell's avatar
gsell committed
389
#endif
390

kraus's avatar
kraus committed
391
    INFOMSG(level3 << "* Computing Preconditioner..." << endl);
392
    IpplTimings::startTimer(FunctionTimer5_m);
frey_m's avatar
frey_m committed
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409
    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;
        }
410
    }
411
    IpplTimings::stopTimer(FunctionTimer5_m);
kraus's avatar
kraus committed
412
    INFOMSG(level3 << "* Done." << endl);
413

gsell's avatar
gsell committed
414 415
    // setup preconditioned iterative solver
    // use old LHS solution as initial guess
kraus's avatar
kraus committed
416
    INFOMSG(level3 << "* Final Setup of Solver..." << endl);
417
    IpplTimings::startTimer(FunctionTimer6_m);
418 419 420
    problem_mp->setOperator(A);
    problem_mp->setLHS(LHS);
    problem_mp->setRHS(RHS);
421 422 423 424 425 426

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

427 428 429
    solver_mp->setProblem( problem_mp);
    if (!problem_mp->isProblemSet()) {
        if (problem_mp->setProblem() == false) {
430
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
431
        }
gsell's avatar
gsell committed
432
    }
433
    IpplTimings::stopTimer(FunctionTimer6_m);
kraus's avatar
kraus committed
434
    INFOMSG(level3 << "* Done." << endl);
435

gsell's avatar
gsell committed
436 437
    double time = MPI_Wtime();

kraus's avatar
kraus committed
438
    INFOMSG(level3 << "* Solving for Space Charge..." << endl);
439
    IpplTimings::startTimer(FunctionTimer7_m);
440
    solver_mp->solve();
frey_m's avatar
frey_m committed
441

442
    IpplTimings::stopTimer(FunctionTimer7_m);
kraus's avatar
kraus committed
443
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
444

kraus's avatar
kraus committed
445 446 447
    std::ofstream timings;
    if (true || verbose_m) {
        time = MPI_Wtime() - time;
448
        double minTime = 0, maxTime = 0, avgTime = 0;
449 450 451
        reduce(time, minTime, 1, std::less<double>());
        reduce(time, maxTime, 1, std::greater<double>());
        reduce(time, avgTime, 1, std::plus<double>());
452 453
        avgTime /= comm_mp->getSize();
        if (comm_mp->getRank() == 0) {
kraus's avatar
kraus committed
454
            char filename[50];
frey_m's avatar
frey_m committed
455 456
            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],
457
                    comm_mp->getSize(), recycleBlocks_m, numBlocks_m, nLHS_m);
frey_m's avatar
frey_m committed
458

kraus's avatar
kraus committed
459
            timings.open(filename, std::ios::app);
460
            timings << solver_mp->getNumIters() << "\t"
kraus's avatar
kraus committed
461 462 463 464 465 466 467 468 469 470 471 472 473 474
                    //<< 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
475
    // Store new LHS in OldLHS
frey_m's avatar
frey_m committed
476 477
    if (nLHS_m > 1) OldLHS.push_front(*(LHS.get()));
    if (OldLHS.size() > nLHS_m) OldLHS.pop_back();
gsell's avatar
gsell committed
478 479

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

frey_m's avatar
frey_m committed
494
    if (itsBunch_m->getLocalTrackStep()+1 == (long long)Track::block->localTimeSteps.front()) {
495 496 497
        A = Teuchos::null;
        LHS = Teuchos::null;
        RHS = Teuchos::null;
498
        prec_mp = Teuchos::null;
499
    }
gsell's avatar
gsell committed
500 501 502
}


503
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
504

gsell's avatar
gsell committed
505
    int NumMyElements = 0;
frey_m's avatar
frey_m committed
506
    std::vector<TpetraGlobalOrdinal_t> MyGlobalElements;
gsell's avatar
gsell committed
507

508
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
509
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
510
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
511 512
                if (bp_m->isInside(idx, idy, idz)) {
                    MyGlobalElements.push_back(bp_m->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
513 514 515
                    NumMyElements++;
                }
            }
kraus's avatar
kraus committed
516
        }
517
    }
frey_m's avatar
frey_m committed
518

519 520 521
    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
522 523
}

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

frey_m's avatar
frey_m committed
526
    A->resumeFill();
527
    A->setAllToScalar(0.0);
528

529 530
    int NumMyElements = map_p->getNodeNumElements();
    auto MyGlobalElements = map_p->getMyGlobalIndices();
gsell's avatar
gsell committed
531

frey_m's avatar
frey_m committed
532 533
    std::vector<TpetraScalar_t> Values(6);
    std::vector<TpetraGlobalOrdinal_t> Indices(6);
gsell's avatar
gsell committed
534

535 536 537
    IrregularDomain::StencilValue_t value;
    IrregularDomain::StencilIndex_t index;

frey_m's avatar
frey_m committed
538
    for (int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
539 540 541

        int NumEntries = 0;

542
        double scaleFactor=1.0;
gsell's avatar
gsell committed
543

544
        bp_m->getBoundaryStencil(MyGlobalElements[i], value, scaleFactor);
545
        RHS->scale(scaleFactor);
gsell's avatar
gsell committed
546

547 548 549 550
        bp_m->getNeighbours(MyGlobalElements[i], index);
        if (index.east != -1) {
            Indices[NumEntries]  = index.east;
            Values[NumEntries++] = value.east;
gsell's avatar
gsell committed
551
        }
552 553 554
        if (index.west != -1) {
            Indices[NumEntries]  = index.west;
            Values[NumEntries++] = value.west;
gsell's avatar
gsell committed
555
        }
556 557 558
        if (index.south != -1) {
            Indices[NumEntries]  = index.south;
            Values[NumEntries++] = value.south;
gsell's avatar
gsell committed
559
        }
560 561 562
        if (index.north != -1) {
            Indices[NumEntries]  = index.north;
            Values[NumEntries++] = value.north;
gsell's avatar
gsell committed
563
        }
564 565 566
        if (index.front != -1) {
            Indices[NumEntries]  = index.front;
            Values[NumEntries++] = value.front;
gsell's avatar
gsell committed
567
        }
568 569 570
        if (index.back != -1) {
            Indices[NumEntries]  = index.back;
            Values[NumEntries++] = value.back;
gsell's avatar
gsell committed
571 572
        }

573
        // if matrix has already been filled (fillComplete()) we can only
gsell's avatar
gsell committed
574
        // replace entries
575

frey_m's avatar
fix bug  
frey_m committed
576
        if (isMatrixfilled_m) {
577
            // off-diagonal entries
578
            A->replaceGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
579
            // diagonal entry
580
            A->replaceGlobalValues(MyGlobalElements[i], 1, &value.center, &MyGlobalElements[i]);
581
        } else {
gsell's avatar
gsell committed
582
            // off-diagonal entries
583
            A->insertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
gsell's avatar
gsell committed
584
            // diagonal entry
585
            A->insertGlobalValues(MyGlobalElements[i], 1, &value.center, &MyGlobalElements[i]);
586
        }
gsell's avatar
gsell committed
587 588
    }

589 590 591
    RCP<ParameterList_t> params = Teuchos::parameterList();
    params->set ("Optimize Storage", true);
    A->fillComplete(params);
frey_m's avatar
fix bug  
frey_m committed
592
    isMatrixfilled_m = true;
gsell's avatar
gsell committed
593 594 595 596 597
}

void MGPoissonSolver::printLoadBalanceStats() {

    //compute some load balance statistics
598 599
    size_t myNumPart = map_p->getNodeNumElements();
    size_t NumPart = map_p->getGlobalNumElements() * 1.0 / comm_mp->getSize();
gsell's avatar
gsell committed
600
    double imbalance = 1.0;
frey_m's avatar
frey_m committed
601
    if (myNumPart >= NumPart)
gsell's avatar
gsell committed
602 603 604 605 606
        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
607 608 609 610 611 612
    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
613

614
    avg /= comm_mp->getSize();
kraus's avatar
kraus committed
615 616
    *gmsg << "LBAL min = " << min << ", max = " << max << ", avg = " << avg << endl;
    *gmsg << "min nr gridpoints = " << minn << ", max nr gridpoints = " << maxn << endl;
gsell's avatar
gsell committed
617 618
}

frey_m's avatar
frey_m committed
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 684 685 686 687 688 689 690 691 692 693
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
694 695 696 697
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;
698
    return os;
gsell's avatar
gsell committed
699
}