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;
210 211
}

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

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

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

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

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

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

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

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

286

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

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

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

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

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

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

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

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

332

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

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

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

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

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

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

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

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

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

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

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

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

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

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


492
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
493

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

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

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

frey_m's avatar
frey_m committed
514
    A->resumeFill();
515
    A->setAllToScalar(0.0);
516

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

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

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

        int NumEntries = 0;

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

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

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

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

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

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

void MGPoissonSolver::printLoadBalanceStats() {

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

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

frey_m's avatar
frey_m committed
605 606 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
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
680 681 682 683
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;
684
    return os;
gsell's avatar
gsell committed
685 686
}

687 688 689 690 691 692 693
// 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: