MGPoissonSolver.cpp 23.9 KB
Newer Older
frey_m's avatar
frey_m committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
//
// Class MGPoissonSolver
//   This class contains methods for solving Poisson's equation for the
//   space charge portion of the calculation.
//
//   A smoothed aggregation based AMG preconditioned iterative solver for space charge
//   \see FFTPoissonSolver
//   \warning This solver is in an EXPERIMENTAL STAGE. For reliable simulations use the FFTPoissonSolver
//
// Copyright (c) 2010 - 2013, Yves Ineichen, ETH Zürich,
//               2013 - 2015, Tülin Kaman, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
frey_m's avatar
frey_m committed
15
// "Toward massively parallel multi-objective optimization with application to
frey_m's avatar
frey_m committed
16 17 18 19 20 21 22 23 24 25 26 27
// particle accelerators" (https://doi.org/10.3929/ethz-a-009792359)
//
// 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/>.
//
28

29
//#define DBG_STENCIL
30

kraus's avatar
kraus committed
31
#include "Solvers/MGPoissonSolver.h"
32 33 34 35 36 37

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

38
#include "Track/Track.h"
gsell's avatar
gsell committed
39 40 41 42 43
#include "Physics/Physics.h"
#include "Utilities/OpalException.h"
#include "Attributes/Attributes.h"
#include "ValueDefinitions/RealVariable.h"
#include "AbstractObjects/OpalData.h"
44
#include "Utilities/Options.h"
45

46 47 48 49 50 51
// #include "Epetra_Operator.h"
// #include "EpetraExt_RowMatrixOut.h"
// #include "Epetra_Import.h"

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

#include "Teuchos_CommandLineProcessor.hpp"

#include "BelosLinearProblem.hpp"
#include "BelosRCGSolMgr.hpp"
57
// #include "BelosEpetraAdapter.hpp"
58 59
#include "BelosBlockCGSolMgr.hpp"

60 61 62
#include <MueLu_CreateTpetraPreconditioner.hpp>
// #include <MueLu_TpetraOperator.hpp>

frey_m's avatar
frey_m committed
63 64 65
// #include "ml_MultiLevelPreconditioner.h"
// #include "ml_MultiLevelOperator.h"
// #include "ml_epetra_utils.h"
66

67 68 69 70
// #include "Isorropia_Exception.hpp"
// #include "Isorropia_Epetra.hpp"
// #include "Isorropia_EpetraRedistributor.hpp"
// #include "Isorropia_EpetraPartitioner.hpp"
71

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

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

78
MGPoissonSolver::MGPoissonSolver ( PartBunch *beam,
79 80 81 82 83 84 85
                                   Mesh_t *mesh,
                                   FieldLayout_t *fl,
                                   std::vector<BoundaryGeometry *> geometries,
                                   std::string itsolver,
                                   std::string interpl,
                                   double tol,
                                   int maxiters,
86
                                   std::string precmode)
frey_m's avatar
fix bug  
frey_m committed
87 88
    : isMatrixfilled_m(false)
    , 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
        hr_m[i] = mesh_m->get_meshSpacing(i);
        orig_nr_m[i] = domain_m[i].length();
    }

frey_m's avatar
frey_m committed
104 105 106
    if (itsolver == "CG") itsolver_m = AZ_cg;
    else if (itsolver == "BICGSTAB") itsolver_m = AZ_bicgstab;
    else if (itsolver == "GMRES") itsolver_m = AZ_gmres;
gsell's avatar
gsell committed
107 108 109
    else throw OpalException("MGPoissonSolver", "No valid iterative solver selected!");

    precmode_m = STD_PREC;
frey_m's avatar
frey_m committed
110 111 112 113
    if (precmode == "STD") precmode_m = STD_PREC;
    else if (precmode == "HIERARCHY") precmode_m = REUSE_HIERARCHY;
    else if (precmode == "REUSE") precmode_m = REUSE_PREC;
    else if (precmode == "NO") precmode_m = NO;
gsell's avatar
gsell committed
114 115 116 117 118 119

    tolerableIterationsCount_m = 2;
    numIter_m = -1;
    forcePreconditionerRecomputation_m = false;

    hasParallelDecompositionChanged_m = true;
120
    repartFreq_m = 1000;
gsell's avatar
gsell committed
121
    useRCB_m = false;
frey_m's avatar
frey_m committed
122
    if (Ippl::Info->getOutputLevel() > 3)
gsell's avatar
gsell committed
123 124 125 126
        verbose_m = true;
    else
        verbose_m = false;

127
    // Find CURRENT geometry
gsell's avatar
gsell committed
128
    currentGeometry = geometries_m[0];
frey_m's avatar
frey_m committed
129 130 131 132 133
    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
134
        } else if (currentGeometry->getTopology() == "BOXCORNER") {
frey_m's avatar
frey_m committed
135 136 137 138 139 140 141 142
            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
143
            bp_m->compute(itsBunch_m->get_hr());
gsell's avatar
gsell committed
144
        } else {
145 146
            throw OpalException("MGPoissonSolver::MGPoissonSolver",
                                "Geometry not known");
gsell's avatar
gsell committed
147
        }
148 149 150 151 152 153 154 155 156 157
    } else {
        NDIndex<3> localId = layout_m->getLocalNDIndex();
        if (localId[0].length() != domain_m[0].length() ||
            localId[1].length() != domain_m[1].length()) {
            throw OpalException("ArbitraryDomain::compute",
                                "The class ArbitraryDomain only works with parallelization\n"
                                "in z-direction.\n"
                                "Please set PARFFTX=FALSE, PARFFTY=FALSE, PARFFTT=TRUE in \n"
                                "the definition of the field solver in the input file.\n");
        }
frey_m's avatar
frey_m committed
158 159
        bp_m = std::unique_ptr<IrregularDomain>(
            new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl));
160
    }
161

162
    map_p = Teuchos::null;
gsell's avatar
gsell committed
163 164 165
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
166 167
    MueLuPrec = Teuchos::null;
    prec_mp = Teuchos::null;
168

gsell's avatar
gsell committed
169 170
    numBlocks_m = Options::numBlocks;
    recycleBlocks_m = Options::recycleBlocks;
171
    nLHS_m = Options::nLHS;
frey_m's avatar
frey_m committed
172
    setupMueLuList();
gsell's avatar
gsell committed
173
    SetupBelosList();
174
    problem_mp = rcp(new Belos::LinearProblem<TpetraScalar_t, TpetraMultiVector_t, TpetraOperator_t>);
175
    // setup Belos solver
frey_m's avatar
frey_m committed
176
    if (numBlocks_m == 0 || recycleBlocks_m == 0)
177
        solver_mp = rcp(new Belos::BlockCGSolMgr<TpetraScalar_t, TpetraMultiVector_t, TpetraOperator_t>());
gsell's avatar
gsell committed
178
    else
179 180 181
        solver_mp = rcp(new Belos::RCGSolMgr<TpetraScalar_t, TpetraMultiVector_t, TpetraOperator_t>());
    solver_mp->setParameters(rcp(&belosList, false));
    convStatusTest = rcp(new Belos::StatusTestGenResNorm<TpetraScalar_t, TpetraMultiVector_t, TpetraOperator_t> (tol));
gsell's avatar
gsell committed
182 183 184
    convStatusTest->defineScaleForm(Belos::NormOfRHS, Belos::TwoNorm);

    //all timers used here
185 186 187 188 189 190 191 192
    FunctionTimer1_m = IpplTimings::getTimer("BGF-IndexCoordMap");
    FunctionTimer2_m = IpplTimings::getTimer("computeMap");
    FunctionTimer3_m = IpplTimings::getTimer("IPPL to RHS");
    FunctionTimer4_m = IpplTimings::getTimer("ComputeStencil");
    FunctionTimer5_m = IpplTimings::getTimer("ML");
    FunctionTimer6_m = IpplTimings::getTimer("Setup");
    FunctionTimer7_m = IpplTimings::getTimer("CG");
    FunctionTimer8_m = IpplTimings::getTimer("LHS to IPPL");
gsell's avatar
gsell committed
193 194
}

195
void MGPoissonSolver::deletePtr() {
196 197
    MueLuPrec = Teuchos::null;
    map_p  = Teuchos::null;
198 199 200
    A      = Teuchos::null;
    LHS    = Teuchos::null;
    RHS    = Teuchos::null;
201
    prec_mp = Teuchos::null;
202 203
}

204 205
MGPoissonSolver::~MGPoissonSolver() {
    deletePtr ();
206 207
    solver_mp = Teuchos::null;
    problem_mp = Teuchos::null;
gsell's avatar
gsell committed
208 209
}

210
void MGPoissonSolver::computePotential(Field_t& /*rho*/, Vector_t /*hr*/, double /*zshift*/) {
gsell's avatar
gsell committed
211 212 213
    throw OpalException("MGPoissonSolver", "not implemented yet");
}

214
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
215
    if (itsBunch_m->getLocalTrackStep()%repartFreq_m == 0){
216
        deletePtr();
frey_m's avatar
frey_m committed
217
        if (useRCB_m)
218
            redistributeWithRCB(localId);
219
        else
220
            IPPLToMap3D(localId);
221 222

        extrapolateLHS();
223 224
    }
}
gsell's avatar
gsell committed
225

226 227 228 229 230
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.
231
    //we also have to redistribute LHS
232

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

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

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

281

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

kraus's avatar
kraus committed
288
    Inform msg("OPAL ", INFORM_ALL_NODES);
gsell's avatar
gsell committed
289 290 291
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
292

frey_m's avatar
frey_m committed
293 294 295
    bp_m->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
    bp_m->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
    bp_m->setNr(nr_m);
296

297 298
    NDIndex<3> localId = layout_m->getLocalNDIndex();

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

    // Define the Map
kraus's avatar
kraus committed
306
    INFOMSG(level3 << "* Computing Map..." << endl);
307
    IpplTimings::startTimer(FunctionTimer2_m);
308
    computeMap(localId);
309
    IpplTimings::stopTimer(FunctionTimer2_m);
kraus's avatar
kraus committed
310
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
311

312
    // Allocate the RHS with the new Epetra Map
313
    if (Teuchos::is_null(RHS))
314 315
        RHS = rcp(new TpetraVector_t(map_p));
    RHS->putScalar(0.0);
316

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

327

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

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

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

gsell's avatar
gsell committed
371 372 373
#ifdef DBG_STENCIL
    EpetraExt::RowMatrixToMatlabFile("DiscrStencil.dat", *A);
#endif
374

kraus's avatar
kraus committed
375
    INFOMSG(level3 << "* Computing Preconditioner..." << endl);
frey_m's avatar
frey_m committed
376
    /* FIXME --> MueLu reuse types
377
    IpplTimings::startTimer(FunctionTimer5_m);
378 379
    if (MueLuPrec == Teuchos::null) {
        MueLuPrec = MueLu::CreateTpetraPreconditioner(A, MLList_m);
380
    } else if (precmode_m == REUSE_HIERARCHY) {
381
        // FIXME MueLuPrec->ReComputePreconditioner();
382
    } else if (precmode_m == REUSE_PREC){
383
    }
384
    IpplTimings::stopTimer(FunctionTimer5_m);
385
    */
kraus's avatar
kraus committed
386
    INFOMSG(level3 << "* Done." << endl);
387

gsell's avatar
gsell committed
388 389
    // setup preconditioned iterative solver
    // use old LHS solution as initial guess
kraus's avatar
kraus committed
390
    INFOMSG(level3 << "* Final Setup of Solver..." << endl);
391
    IpplTimings::startTimer(FunctionTimer6_m);
392 393 394
    problem_mp->setOperator(A);
    problem_mp->setLHS(LHS);
    problem_mp->setRHS(RHS);
frey_m's avatar
frey_m committed
395 396 397 398
    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);
    }
399 400 401 402
    problem_mp->setLeftPrec(prec_mp);
    solver_mp->setProblem( problem_mp);
    if (!problem_mp->isProblemSet()) {
        if (problem_mp->setProblem() == false) {
403
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
404
        }
gsell's avatar
gsell committed
405
    }
406
    IpplTimings::stopTimer(FunctionTimer6_m);
kraus's avatar
kraus committed
407
    INFOMSG(level3 << "* Done." << endl);
408

gsell's avatar
gsell committed
409 410
    double time = MPI_Wtime();

kraus's avatar
kraus committed
411
    INFOMSG(level3 << "* Solving for Space Charge..." << endl);
412
    IpplTimings::startTimer(FunctionTimer7_m);
413
    solver_mp->solve();
414
    IpplTimings::stopTimer(FunctionTimer7_m);
kraus's avatar
kraus committed
415
    INFOMSG(level3 << "* Done." << endl);
gsell's avatar
gsell committed
416

kraus's avatar
kraus committed
417 418 419
    std::ofstream timings;
    if (true || verbose_m) {
        time = MPI_Wtime() - time;
420
        double minTime = 0, maxTime = 0, avgTime = 0;
421 422 423
        reduce(time, minTime, 1, std::less<double>());
        reduce(time, maxTime, 1, std::greater<double>());
        reduce(time, avgTime, 1, std::plus<double>());
424 425
        avgTime /= comm_mp->getSize();
        if (comm_mp->getRank() == 0) {
kraus's avatar
kraus committed
426
            char filename[50];
frey_m's avatar
frey_m committed
427 428
            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],
429
                    comm_mp->getSize(), recycleBlocks_m, numBlocks_m, nLHS_m);
frey_m's avatar
frey_m committed
430

kraus's avatar
kraus committed
431
            timings.open(filename, std::ios::app);
432
            timings << solver_mp->getNumIters() << "\t"
kraus's avatar
kraus committed
433 434 435 436 437 438 439 440 441 442 443 444 445 446
                    //<< 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
447
    // Store new LHS in OldLHS
frey_m's avatar
frey_m committed
448 449
    if (nLHS_m > 1) OldLHS.push_front(*(LHS.get()));
    if (OldLHS.size() > nLHS_m) OldLHS.pop_back();
gsell's avatar
gsell committed
450 451

    //now transfer solution back to IPPL grid
452 453
    IpplTimings::startTimer(FunctionTimer8_m);
    id = 0;
gsell's avatar
gsell committed
454
    rho = 0.0;
455
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
456
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
kraus's avatar
kraus committed
457
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
458
                  NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
frey_m's avatar
frey_m committed
459
                  if (bp_m->isInside(idx, idy, idz))
460
                     rho.localElement(l) = LHS->getData()[id++] * scaleFactor;
gsell's avatar
gsell committed
461 462 463
            }
        }
    }
464
    IpplTimings::stopTimer(FunctionTimer8_m);
465

frey_m's avatar
frey_m committed
466
    if (itsBunch_m->getLocalTrackStep()+1 == (long long)Track::block->localTimeSteps.front()) {
467 468 469
        A = Teuchos::null;
        LHS = Teuchos::null;
        RHS = Teuchos::null;
470
        prec_mp = Teuchos::null;
471
    }
gsell's avatar
gsell committed
472 473 474
}


475
void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
gsell's avatar
gsell committed
476 477 478

    int numMyGridPoints = 0;

479
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
480
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
kraus's avatar
kraus committed
481
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
482
                if (bp_m->isInside(idx, idy, idz))
gsell's avatar
gsell committed
483
                    numMyGridPoints++;
484
             }
gsell's avatar
gsell committed
485
        }
486
     }
gsell's avatar
gsell committed
487 488


489 490 491 492 493
    int indexbase = 0;
    Teuchos::RCP<TpetraMap_t> bmap = Teuchos::rcp(new TpetraMap_t(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
                                                                  numMyGridPoints, indexbase, comm_mp));
    Teuchos::RCP<const TpetraMultiVector_t> coords = Teuchos::rcp(
        new TpetraMultiVector_t(bmap, 3, false));
gsell's avatar
gsell committed
494

frey_m's avatar
frey_m committed
495
/*FIXME --> Compare src/Solvers/AMR_MG/MueLuPreconditioner.hpp
496 497 498
    double *v;
    int stride = 0, stride2 = 0;
    coords->extractView(&v, &stride);
gsell's avatar
gsell committed
499 500
    stride2 = 2 * stride;

501
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
502
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
kraus's avatar
kraus committed
503
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
504
                if (bp_m->isInside(idx, idy, idz)) {
505 506 507
                    v[0] = (double)idx;
                    v[stride] = (double)idy;
                    v[stride2] = (double)idz;
gsell's avatar
gsell committed
508 509 510 511 512
                    v++;
                }
            }
        }
    }
513
    */
gsell's avatar
gsell committed
514 515 516 517 518 519 520

    Teuchos::ParameterList paramlist;
    paramlist.set("Partitioning Method", "RCB");
    Teuchos::ParameterList &sublist = paramlist.sublist("ZOLTAN");
    sublist.set("RCB_RECTILINEAR_BLOCKS", "1");
    sublist.set("DEBUG_LEVEL", "1");

521 522
    /*
     * FIXME
frey_m's avatar
frey_m committed
523 524 525
    Teuchos::RCP<Isorropia::Epetra::Partitioner> part = Teuchos::rcp(
        new Isorropia::Epetra::Partitioner(coords, paramlist));

gsell's avatar
gsell committed
526
    Isorropia::Epetra::Redistributor rd(part);
527
    Teuchos::RCP<TpetraMultiVector_t> newcoords = rd.redistribute(*coords);
gsell's avatar
gsell committed
528

529
    newcoords->extractView(&v, &stride);
gsell's avatar
gsell committed
530 531
    stride2 = 2 * stride;
    numMyGridPoints = 0;
532
    */
gsell's avatar
gsell committed
533
    std::vector<int> MyGlobalElements;
534 535
    /* FIXME
    for (unsigned int i = 0; i < newcoords->getLocalLength(); i++) {
frey_m's avatar
frey_m committed
536
        MyGlobalElements.push_back(bp_m->getIdx(v[0], v[stride], v[stride2]));
gsell's avatar
gsell committed
537 538 539
        v++;
        numMyGridPoints++;
    }
540 541 542
*/
    map_p = Teuchos::rcp(new TpetraMap_t(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
                                         &MyGlobalElements[0], numMyGridPoints, indexbase, comm_mp));
gsell's avatar
gsell committed
543 544
}

545
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
546

gsell's avatar
gsell committed
547
    int NumMyElements = 0;
548
    std::vector<int> MyGlobalElements;
gsell's avatar
gsell committed
549

550
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
551
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
552
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
frey_m's avatar
frey_m committed
553 554
                if (bp_m->isInside(idx, idy, idz)) {
                    MyGlobalElements.push_back(bp_m->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
555 556 557
                    NumMyElements++;
                }
            }
kraus's avatar
kraus committed
558
        }
559
    }
560 561 562
    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
563 564
}

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

frey_m's avatar
frey_m committed
567
    A->resumeFill();
568
    A->setAllToScalar(0.0);
569

570 571
    int NumMyElements = map_p->getNodeNumElements();
    auto MyGlobalElements = map_p->getMyGlobalIndices();
gsell's avatar
gsell committed
572

573 574
    std::vector<double> Values(6);
    std::vector<int> Indices(6);
gsell's avatar
gsell committed
575

frey_m's avatar
frey_m committed
576
    for (int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
577 578 579

        int NumEntries = 0;

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

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

frey_m's avatar
frey_m committed
586
        bp_m->getNeighbours(MyGlobalElements[i], W, E, S, N, F, B);
frey_m's avatar
frey_m committed
587
        if (E != -1) {
gsell's avatar
gsell committed
588 589 590
            Indices[NumEntries] = E;
            Values[NumEntries++] = EV;
        }
frey_m's avatar
frey_m committed
591
        if (W != -1) {
gsell's avatar
gsell committed
592 593 594
            Indices[NumEntries] = W;
            Values[NumEntries++] = WV;
        }
frey_m's avatar
frey_m committed
595
        if (S != -1) {
gsell's avatar
gsell committed
596 597 598
            Indices[NumEntries] = S;
            Values[NumEntries++] = SV;
        }
frey_m's avatar
frey_m committed
599
        if (N != -1) {
gsell's avatar
gsell committed
600 601 602
            Indices[NumEntries] = N;
            Values[NumEntries++] = NV;
        }
frey_m's avatar
frey_m committed
603
        if (F != -1) {
gsell's avatar
gsell committed
604 605 606
            Indices[NumEntries] = F;
            Values[NumEntries++] = FV;
        }
frey_m's avatar
frey_m committed
607
        if (B != -1) {
gsell's avatar
gsell committed
608 609 610 611
            Indices[NumEntries] = B;
            Values[NumEntries++] = BV;
        }

612
        // if matrix has already been filled (fillComplete()) we can only
gsell's avatar
gsell committed
613
        // replace entries
614

frey_m's avatar
fix bug  
frey_m committed
615
        if (isMatrixfilled_m) {
616
            // off-diagonal entries
617
            A->replaceGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
618
            // diagonal entry
619
            A->replaceGlobalValues(MyGlobalElements[i], 1, &CV, &MyGlobalElements[i]);
620
        } else {
gsell's avatar
gsell committed
621
            // off-diagonal entries
622
            A->insertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
gsell's avatar
gsell committed
623
            // diagonal entry
624
            A->insertGlobalValues(MyGlobalElements[i], 1, &CV, &MyGlobalElements[i]);
625
        }
gsell's avatar
gsell committed
626 627
    }

628 629 630
    RCP<ParameterList_t> params = Teuchos::parameterList();
    params->set ("Optimize Storage", true);
    A->fillComplete(params);
frey_m's avatar
fix bug  
frey_m committed
631
    isMatrixfilled_m = true;
gsell's avatar
gsell committed
632 633 634 635 636
}

void MGPoissonSolver::printLoadBalanceStats() {

    //compute some load balance statistics
637 638
    size_t myNumPart = map_p->getNodeNumElements();
    size_t NumPart = map_p->getGlobalNumElements() * 1.0 / comm_mp->getSize();
gsell's avatar
gsell committed
639
    double imbalance = 1.0;
frey_m's avatar
frey_m committed
640
    if (myNumPart >= NumPart)
gsell's avatar
gsell committed
641 642 643 644 645 646 647 648 649 650 651 652
        imbalance += (myNumPart - NumPart) / NumPart;
    else
        imbalance += (NumPart - myNumPart) / NumPart;

    double max = 0.0, min = 0.0, avg = 0.0;
    int minn = 0, maxn = 0;
    MPI_Reduce(&imbalance, &min, 1, MPI_DOUBLE, MPI_MIN, 0, Ippl::getComm());
    MPI_Reduce(&imbalance, &max, 1, MPI_DOUBLE, MPI_MAX, 0, Ippl::getComm());
    MPI_Reduce(&imbalance, &avg, 1, MPI_DOUBLE, MPI_SUM, 0, Ippl::getComm());
    MPI_Reduce(&myNumPart, &minn, 1, MPI_INT, MPI_MIN, 0, Ippl::getComm());
    MPI_Reduce(&myNumPart, &maxn, 1, MPI_INT, MPI_MAX, 0, Ippl::getComm());

653
    avg /= comm_mp->getSize();
kraus's avatar
kraus committed
654 655
    *gmsg << "LBAL min = " << min << ", max = " << max << ", avg = " << avg << endl;
    *gmsg << "min nr gridpoints = " << minn << ", max nr gridpoints = " << maxn << endl;
gsell's avatar
gsell committed
656 657 658 659 660 661
}

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;
662
    return os;
gsell's avatar
gsell committed
663 664
}

665 666 667 668 669 670 671
// 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: