MGPoissonSolver.cpp 25.4 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9
//
// Class MGPoissonSolver
//   This class contains methods for solving Poisson's equation for the
//   space charge portion of the calculation.
//
//   A smoothed aggregation based AMG preconditioned iterative solver for space charge
//   \see FFTPoissonSolver
//   \warning This solver is in an EXPERIMENTAL STAGE. For reliable simulations use the FFTPoissonSolver
//
10 11 12
// Copyright (c) 2008,        Yves Ineichen, ETH Zürich,
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
//               2017 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
frey_m's avatar
frey_m committed
13 14
// All rights reserved
//
15 16 17 18 19
// Implemented as part of the master thesis
// "A Parallel Multigrid Solver for Beam Dynamics"
// and the paper
// "A fast parallel Poisson solver on irregular domains applied to beam dynamics simulations"
// (https://doi.org/10.1016/j.jcp.2010.02.022)
frey_m's avatar
frey_m committed
20
//
frey_m's avatar
frey_m committed
21
// In 2020, the code was ported to the second generation Trilinos packages,
frey_m's avatar
frey_m committed
22
// i.e., Epetra --> Tpetra, ML --> MueLu. See also issue #507.
frey_m's avatar
frey_m committed
23
//
frey_m's avatar
frey_m committed
24 25 26 27 28 29 30 31 32 33
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
34

35
//#define DBG_STENCIL
36

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

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

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

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

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

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

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

69
#include <MueLu_CreateTpetraPreconditioner.hpp>
70

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

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

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

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

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

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

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

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

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

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

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

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

gsell's avatar
gsell committed
200 201

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

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

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

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

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

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

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

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

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

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

295

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

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

frey_m's avatar
frey_m committed
307 308 309
    bp_m->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
    bp_m->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
    bp_m->setNr(nr_m);
310

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

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

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

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

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

341

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

        int NumEntries = 0;

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

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

frey_m's avatar
frey_m committed
546
        bp_m->getNeighbours(MyGlobalElements[i], W, E, S, N, F, B);
frey_m's avatar
frey_m committed
547
        if (E != -1) {
gsell's avatar
gsell committed
548 549 550
            Indices[NumEntries] = E;
            Values[NumEntries++] = EV;
        }
frey_m's avatar
frey_m committed
551
        if (W != -1) {
gsell's avatar
gsell committed
552 553 554
            Indices[NumEntries] = W;
            Values[NumEntries++] = WV;
        }
frey_m's avatar
frey_m committed
555
        if (S != -1) {
gsell's avatar
gsell committed
556 557 558
            Indices[NumEntries] = S;
            Values[NumEntries++] = SV;
        }
frey_m's avatar
frey_m committed
559
        if (N != -1) {
gsell's avatar
gsell committed
560 561 562
            Indices[NumEntries] = N;
            Values[NumEntries++] = NV;
        }
frey_m's avatar
frey_m committed
563
        if (F != -1) {
gsell's avatar
gsell committed
564 565 566
            Indices[NumEntries] = F;
            Values[NumEntries++] = FV;
        }
frey_m's avatar
frey_m committed
567
        if (B != -1) {
gsell's avatar
gsell committed
568 569 570 571
            Indices[NumEntries] = B;
            Values[NumEntries++] = BV;
        }

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

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

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

void MGPoissonSolver::printLoadBalanceStats() {

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

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

frey_m's avatar
frey_m committed
618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692
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
693 694 695 696
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;
697
    return os;
gsell's avatar
gsell committed
698 699
}

700 701 702 703 704 705 706
// 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: