MGPoissonSolver.cpp 25.1 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));
frey_m's avatar
frey_m committed
129
            bp_m->compute(itsBunch_m->get_hr());
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
        bp_m = std::unique_ptr<IrregularDomain>(
            new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl));
146
    }
147

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

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

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

gsell's avatar
gsell committed
194 195

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

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

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

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

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

        extrapolateLHS();
231 232
    }
}
gsell's avatar
gsell committed
233

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

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

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

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

289

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

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

frey_m's avatar
frey_m committed
301 302 303
    bp_m->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
    bp_m->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
    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

frey_m's avatar
frey_m committed
530
    for (int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
531 532 533

        int NumEntries = 0;

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

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

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

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

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

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

void MGPoissonSolver::printLoadBalanceStats() {

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

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

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

694 695 696 697 698 699 700
// 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: