MGPoissonSolver.cpp 24.8 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 84
    : isMatrixfilled_m(false)
    , geometries_m(geometries)
85 86 87 88 89 90 91
    , tol_m(tol)
    , maxiters_m(maxiters)
    , comm_mp(new Comm_t(Ippl::getComm()))
    , itsBunch_m(beam)
    , mesh_m(mesh)
    , layout_m(fl)
{
92

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

95
    for (int i = 0; i < 3; i++) {
gsell's avatar
gsell committed
96 97 98 99 100
        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
101
    if (precmode == "HIERARCHY") precmode_m = REUSE_HIERARCHY;
frey_m's avatar
frey_m committed
102 103
    else if (precmode == "REUSE") precmode_m = REUSE_PREC;
    else if (precmode == "NO") precmode_m = NO;
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 136 137 138 139 140 141
    } else {
        NDIndex<3> localId = layout_m->getLocalNDIndex();
        if (localId[0].length() != domain_m[0].length() ||
            localId[1].length() != domain_m[1].length()) {
            throw OpalException("ArbitraryDomain::compute",
                                "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 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188
    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") {
        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!");
    }

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

gsell's avatar
gsell committed
191 192

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

203
void MGPoissonSolver::deletePtr() {
204
    map_p  = Teuchos::null;
205 206 207
    A      = Teuchos::null;
    LHS    = Teuchos::null;
    RHS    = Teuchos::null;
208
    prec_mp = Teuchos::null;
209 210
}

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

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

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

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

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

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

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

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

285

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

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

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

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

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

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

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

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

331

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

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

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

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

kraus's avatar
kraus committed
380
    INFOMSG(level3 << "* Computing Preconditioner..." << endl);
381
    IpplTimings::startTimer(FunctionTimer5_m);
frey_m's avatar
frey_m committed
382 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 NO:
        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 411 412 413 414
    problem_mp->setOperator(A);
    problem_mp->setLHS(LHS);
    problem_mp->setRHS(RHS);
    problem_mp->setLeftPrec(prec_mp);
    solver_mp->setProblem( problem_mp);
    if (!problem_mp->isProblemSet()) {
        if (problem_mp->setProblem() == false) {
415
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
416
        }
gsell's avatar
gsell committed
417
    }
418
    IpplTimings::stopTimer(FunctionTimer6_m);
kraus's avatar
kraus committed
419
    INFOMSG(level3 << "* Done." << endl);
420

gsell's avatar
gsell committed
421 422
    double time = MPI_Wtime();

kraus's avatar
kraus committed
423
    INFOMSG(level3 << "* Solving for Space Charge..." << endl);
424
    IpplTimings::startTimer(FunctionTimer7_m);
425
    solver_mp->solve();
426
    IpplTimings::stopTimer(FunctionTimer7_m);
kraus's avatar
kraus committed
427
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
428

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

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

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

frey_m's avatar
frey_m committed
478
    if (itsBunch_m->getLocalTrackStep()+1 == (long long)Track::block->localTimeSteps.front()) {
479 480 481
        A = Teuchos::null;
        LHS = Teuchos::null;
        RHS = Teuchos::null;
482
        prec_mp = Teuchos::null;
483
    }
gsell's avatar
gsell committed
484 485 486
}


487
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
488

gsell's avatar
gsell committed
489
    int NumMyElements = 0;
490
    std::vector<int> MyGlobalElements;
gsell's avatar
gsell committed
491

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

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

frey_m's avatar
frey_m committed
509
    A->resumeFill();
510
    A->setAllToScalar(0.0);
511

512 513
    int NumMyElements = map_p->getNodeNumElements();
    auto MyGlobalElements = map_p->getMyGlobalIndices();
gsell's avatar
gsell committed
514

515 516
    std::vector<double> Values(6);
    std::vector<int> Indices(6);
gsell's avatar
gsell committed
517

frey_m's avatar
frey_m committed
518
    for (int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
519 520 521

        int NumEntries = 0;

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

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

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

554
        // if matrix has already been filled (fillComplete()) we can only
gsell's avatar
gsell committed
555
        // replace entries
556

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

570 571 572
    RCP<ParameterList_t> params = Teuchos::parameterList();
    params->set ("Optimize Storage", true);
    A->fillComplete(params);
frey_m's avatar
fix bug  
frey_m committed
573
    isMatrixfilled_m = true;
gsell's avatar
gsell committed
574 575 576 577 578
}

void MGPoissonSolver::printLoadBalanceStats() {

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

595
    avg /= comm_mp->getSize();
kraus's avatar
kraus committed
596 597
    *gmsg << "LBAL min = " << min << ", max = " << max << ", avg = " << avg << endl;
    *gmsg << "min nr gridpoints = " << minn << ", max nr gridpoints = " << maxn << endl;
gsell's avatar
gsell committed
598 599
}

frey_m's avatar
frey_m committed
600 601 602 603 604 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
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 NO:
        case STD_PREC:
        default:
            MueLuList_m.set("reuse: type", "none");
            break;
    }
}

gsell's avatar
gsell committed
676 677 678 679
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;
680
    return os;
gsell's avatar
gsell committed
681 682
}

683 684 685 686 687 688 689
// 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: