MGPoissonSolver.cpp 25.3 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 43

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

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

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

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

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

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

68
#include <MueLu_CreateTpetraPreconditioner.hpp>
69

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

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

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

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

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

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

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

150
    map_p = Teuchos::null;
gsell's avatar
gsell committed
151 152 153
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
154
    prec_mp = Teuchos::null;
155

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

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

gsell's avatar
gsell committed
196 197

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

208
void MGPoissonSolver::deletePtr() {
209
    map_p  = Teuchos::null;
210 211 212
    A      = Teuchos::null;
    LHS    = Teuchos::null;
    RHS    = Teuchos::null;
213
    prec_mp = Teuchos::null;
frey_m's avatar
frey_m committed
214
    isMatrixfilled_m = false;
215 216
}

217 218
MGPoissonSolver::~MGPoissonSolver() {
    deletePtr ();
219 220
    solver_mp = Teuchos::null;
    problem_mp = Teuchos::null;
gsell's avatar
gsell committed
221 222
}

223
void MGPoissonSolver::computePotential(Field_t& /*rho*/, Vector_t /*hr*/, double /*zshift*/) {
gsell's avatar
gsell committed
224 225 226
    throw OpalException("MGPoissonSolver", "not implemented yet");
}

227
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
228
    if (itsBunch_m->getLocalTrackStep()%repartFreq_m == 0){
229
        deletePtr();
frey_m's avatar
frey_m committed
230
        IPPLToMap3D(localId);
231 232

        extrapolateLHS();
233 234
    }
}
gsell's avatar
gsell committed
235

236 237 238 239 240
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.
241
    //we also have to redistribute LHS
242

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

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

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

291

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

kraus's avatar
kraus committed
298
    Inform msg("OPAL ", INFORM_ALL_NODES);
gsell's avatar
gsell committed
299 300 301
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
302

frey_m's avatar
frey_m committed
303
    bp_m->setNr(nr_m);
304

305 306
    NDIndex<3> localId = layout_m->getLocalNDIndex();

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

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

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

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

335

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

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

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

gsell's avatar
gsell committed
381
#ifdef DBG_STENCIL
frey_m's avatar
frey_m committed
382 383
    Tpetra::MatrixMarket::Writer<TpetraCrsMatrix_t>::writeSparseFile(
        "DiscrStencil.dat", A);
gsell's avatar
gsell committed
384
#endif
385

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

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

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

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

gsell's avatar
gsell committed
431 432
    double time = MPI_Wtime();

kraus's avatar
kraus committed
433
    INFOMSG(level3 << "* Solving for Space Charge..." << endl);
434
    IpplTimings::startTimer(FunctionTimer7_m);
435
    solver_mp->solve();
frey_m's avatar
frey_m committed
436

437
    IpplTimings::stopTimer(FunctionTimer7_m);
kraus's avatar
kraus committed
438
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
439

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

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

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

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


498
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
499

gsell's avatar
gsell committed
500
    int NumMyElements = 0;
frey_m's avatar
frey_m committed
501
    std::vector<TpetraGlobalOrdinal_t> MyGlobalElements;
gsell's avatar
gsell committed
502

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

514 515 516
    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
517 518
}

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

frey_m's avatar
frey_m committed
521
    A->resumeFill();
522
    A->setAllToScalar(0.0);
523

524 525
    int NumMyElements = map_p->getNodeNumElements();
    auto MyGlobalElements = map_p->getMyGlobalIndices();
gsell's avatar
gsell committed
526

frey_m's avatar
frey_m committed
527 528
    std::vector<TpetraScalar_t> Values(6);
    std::vector<TpetraGlobalOrdinal_t> Indices(6);
gsell's avatar
gsell committed
529

530 531 532
    IrregularDomain::StencilValue_t value;
    IrregularDomain::StencilIndex_t index;

frey_m's avatar
frey_m committed
533
    for (int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
534 535 536

        int NumEntries = 0;

537
        double scaleFactor=1.0;
gsell's avatar
gsell committed
538

539
        bp_m->getBoundaryStencil(MyGlobalElements[i], value, scaleFactor);
540
        RHS->scale(scaleFactor);
gsell's avatar
gsell committed
541

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

568
        // if matrix has already been filled (fillComplete()) we can only
gsell's avatar
gsell committed
569
        // replace entries
570

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

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

void MGPoissonSolver::printLoadBalanceStats() {

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

609
    avg /= comm_mp->getSize();
kraus's avatar
kraus committed
610 611
    *gmsg << "LBAL min = " << min << ", max = " << max << ", avg = " << avg << endl;
    *gmsg << "min nr gridpoints = " << minn << ", max nr gridpoints = " << maxn << endl;
gsell's avatar
gsell committed
612 613
}

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

696 697 698 699 700 701 702
// 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: