MGPoissonSolver.cpp 24.9 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10
//
// 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
//
// Copyright (c) 2010 - 2013, Yves Ineichen, ETH Zürich,
frey_m's avatar
frey_m committed
11 12
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland,
//                      2020, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
13 14 15
// All rights reserved
//
// Implemented as part of the PhD thesis
frey_m's avatar
frey_m committed
16
// "Toward massively parallel multi-objective optimization with application to
frey_m's avatar
frey_m committed
17 18
// particle accelerators" (https://doi.org/10.3929/ethz-a-009792359)
//
frey_m's avatar
frey_m committed
19
// In 2020, the code was ported to the second generation Trilinos packages,
frey_m's avatar
frey_m committed
20
// i.e., Epetra --> Tpetra, ML --> MueLu. See also issue #507.
frey_m's avatar
frey_m committed
21
//
frey_m's avatar
frey_m committed
22 23 24 25 26 27 28 29 30 31
// 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/>.
//
32

33
//#define DBG_STENCIL
34

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

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

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

50 51
#include <Tpetra_Import.hpp>
#include <BelosTpetraAdapter.hpp>
52

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

57 58 59 60 61
#include "Teuchos_CommandLineProcessor.hpp"

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

66
#include <MueLu_CreateTpetraPreconditioner.hpp>
67

gsell's avatar
gsell committed
68 69
#include <algorithm>

70 71
using Teuchos::RCP;
using Teuchos::rcp;
gsell's avatar
gsell committed
72 73
using Physics::c;

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

gsell's avatar
gsell committed
94 95
    domain_m = layout_m->getDomain();

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

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

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

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

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

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

gsell's avatar
gsell committed
192 193

    //all timers used here
194 195 196 197
    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
198
    FunctionTimer5_m = IpplTimings::getTimer("MueLu");
199 200 201
    FunctionTimer6_m = IpplTimings::getTimer("Setup");
    FunctionTimer7_m = IpplTimings::getTimer("CG");
    FunctionTimer8_m = IpplTimings::getTimer("LHS to IPPL");
gsell's avatar
gsell committed
202 203
}

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

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

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

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

        extrapolateLHS();
229 230
    }
}
gsell's avatar
gsell committed
231

232 233 234 235 236
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.
237
    //we also have to redistribute LHS
238

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

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

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

287

gsell's avatar
gsell committed
288 289
// 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
290
// XXX: use matrix stencil in computation directly (no Tpetra, define operators
gsell's avatar
gsell committed
291 292
// on IPPL GRID)
void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr) {
293

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

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

303 304
    NDIndex<3> localId = layout_m->getLocalNDIndex();

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

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

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

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

333

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

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

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

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

kraus's avatar
kraus committed
382
    INFOMSG(level3 << "* Computing Preconditioner..." << endl);
383
    IpplTimings::startTimer(FunctionTimer5_m);
frey_m's avatar
frey_m committed
384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400
    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;
        }
401
    }
402
    IpplTimings::stopTimer(FunctionTimer5_m);
kraus's avatar
kraus committed
403
    INFOMSG(level3 << "* Done." << endl);
404

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

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

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

gsell's avatar
gsell committed
427 428
    double time = MPI_Wtime();

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

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

kraus's avatar
kraus committed
436 437 438
    std::ofstream timings;
    if (true || verbose_m) {
        time = MPI_Wtime() - time;
439
        double minTime = 0, maxTime = 0, avgTime = 0;
440 441 442
        reduce(time, minTime, 1, std::less<double>());
        reduce(time, maxTime, 1, std::greater<double>());
        reduce(time, avgTime, 1, std::plus<double>());
443 444
        avgTime /= comm_mp->getSize();
        if (comm_mp->getRank() == 0) {
kraus's avatar
kraus committed
445
            char filename[50];
frey_m's avatar
frey_m committed
446 447
            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],
448
                    comm_mp->getSize(), recycleBlocks_m, numBlocks_m, nLHS_m);
frey_m's avatar
frey_m committed
449

kraus's avatar
kraus committed
450
            timings.open(filename, std::ios::app);
451
            timings << solver_mp->getNumIters() << "\t"
kraus's avatar
kraus committed
452 453 454 455 456 457 458 459 460 461 462 463 464 465
                    //<< 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
466
    // Store new LHS in OldLHS
frey_m's avatar
frey_m committed
467 468
    if (nLHS_m > 1) OldLHS.push_front(*(LHS.get()));
    if (OldLHS.size() > nLHS_m) OldLHS.pop_back();
gsell's avatar
gsell committed
469 470

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

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


494
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
495

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

499
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
500
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
501
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
502 503
                if (bp_m->isInside(idx, idy, idz)) {
                    MyGlobalElements.push_back(bp_m->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
504 505 506
                    NumMyElements++;
                }
            }
kraus's avatar
kraus committed
507
        }
508
    }
509 510 511
    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
512 513
}

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

frey_m's avatar
frey_m committed
516
    A->resumeFill();
517
    A->setAllToScalar(0.0);
518

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

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

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

        int NumEntries = 0;

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

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

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

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

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

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

void MGPoissonSolver::printLoadBalanceStats() {

    //compute some load balance statistics
586 587
    size_t myNumPart = map_p->getNodeNumElements();
    size_t NumPart = map_p->getGlobalNumElements() * 1.0 / comm_mp->getSize();
gsell's avatar
gsell committed
588
    double imbalance = 1.0;
frey_m's avatar
frey_m committed
589
    if (myNumPart >= NumPart)
gsell's avatar
gsell committed
590 591 592 593 594
        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
595 596 597 598 599 600
    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
601

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

frey_m's avatar
frey_m committed
607 608 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
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
682 683 684 685
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;
686
    return os;
gsell's avatar
gsell committed
687 688
}

689 690 691 692 693 694 695
// 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: