MGPoissonSolver.cpp 19.9 KB
Newer Older
1
#ifdef HAVE_SAAMG_SOLVER
gsell's avatar
gsell committed
2
#define DBG_STENCIL
3

kraus's avatar
kraus committed
4
#include "Solvers/MGPoissonSolver.h"
5 6 7 8 9 10 11 12

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

#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
13 14 15 16 17 18 19 20
#include "Physics/Physics.h"
#include "Attributes/Attributes.h"
#include "ValueDefinitions/RealVariable.h"
#include "AbstractObjects/OpalData.h"
#include "Utilities/OpalException.h"
#include "Attributes/Attributes.h"
#include "ValueDefinitions/RealVariable.h"
#include "AbstractObjects/OpalData.h"
kraus's avatar
kraus committed
21
#include "Utilities/OpalOptions.h"
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
#include "Epetra_Map.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Operator.h"
#include "EpetraExt_RowMatrixOut.h"
#include "Epetra_Import.h"

#include "Teuchos_CommandLineProcessor.hpp"

#include "BelosLinearProblem.hpp"
#include "BelosRCGSolMgr.hpp"
#include "BelosEpetraAdapter.hpp"
#include "BelosBlockCGSolMgr.hpp"

#include "ml_MultiLevelPreconditioner.h"
#include "ml_MultiLevelOperator.h"
#include "ml_epetra_utils.h"

#include "Isorropia_Exception.hpp"
#include "Isorropia_Epetra.hpp"
#include "Isorropia_EpetraRedistributor.hpp"
#include "Isorropia_EpetraPartitioner.hpp"

#pragma GCC diagnostic pop

gsell's avatar
gsell committed
50
#include <algorithm>
51
#include <omp.h>
gsell's avatar
gsell committed
52

53 54
using Teuchos::RCP;
using Teuchos::rcp;
gsell's avatar
gsell committed
55 56
using Physics::c;

57
MGPoissonSolver::MGPoissonSolver ( PartBunch &beam,
58 59 60 61 62 63 64
                                   Mesh_t *mesh,
                                   FieldLayout_t *fl,
                                   std::vector<BoundaryGeometry *> geometries,
                                   std::string itsolver,
                                   std::string interpl,
                                   double tol,
                                   int maxiters,
65
                                   std::string precmode):
66
    itsBunch_m(&beam),
67 68 69 70 71
    mesh_m(mesh),
    layout_m(fl),
    geometries_m(geometries),
    tol_m(tol),
    maxiters_m(maxiters),
72
    Comm(Ippl::getComm()) {
73

gsell's avatar
gsell committed
74 75
    domain_m = layout_m->getDomain();

76
    for (int i = 0; i < 3; i++) {
gsell's avatar
gsell committed
77 78 79 80 81 82 83 84 85 86 87 88 89
        hr_m[i] = mesh_m->get_meshSpacing(i);
        orig_nr_m[i] = domain_m[i].length();
    }

    if(itsolver == "CG") itsolver_m = AZ_cg;
    else if(itsolver == "BiCGStab") itsolver_m = AZ_bicgstab;
    else if(itsolver == "GMRES") itsolver_m = AZ_gmres;
    else throw OpalException("MGPoissonSolver", "No valid iterative solver selected!");

    precmode_m = STD_PREC;
    if(precmode == "STD") precmode_m = STD_PREC;
    else if(precmode == "HIERARCHY") precmode_m = REUSE_HIERARCHY;
    else if(precmode == "REUSE") precmode_m = REUSE_PREC;
90
    else if(precmode == "NO") precmode_m = NO;
gsell's avatar
gsell committed
91 92 93 94 95 96 97

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

    hasParallelDecompositionChanged_m = true;
    useRCB_m = false;
98
    if(Ippl::Info->getOutputLevel() > 1)
gsell's avatar
gsell committed
99 100 101 102
        verbose_m = true;
    else
        verbose_m = false;

103
    // Find CURRENT geometry
gsell's avatar
gsell committed
104 105
    currentGeometry = geometries_m[0];
    if(currentGeometry->getFilename() == "") {
106
        if(currentGeometry->getTopology() == "ELLIPTIC"){
gsell's avatar
gsell committed
107
            bp = new EllipticDomain(currentGeometry->getA(), currentGeometry->getB(), orig_nr_m, hr_m, interpl);
108
	} else if (currentGeometry->getTopology() == "BOXCORNER") {
109
            bp = new BoxCornerDomain(currentGeometry->getA(), currentGeometry->getB(), currentGeometry->getC(), currentGeometry->getLength(),currentGeometry->getL1(), currentGeometry->getL2(), orig_nr_m, hr_m, interpl);
110
            bp->Compute(itsBunch_m->get_hr());
gsell's avatar
gsell committed
111
        } else {
112 113
            throw OpalException("MGPoissonSolver::MGPoissonSolver",
                                "Geometry not known");
gsell's avatar
gsell committed
114
        }
115
    } else
116 117
	bp = new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl);

gsell's avatar
gsell committed
118 119 120 121
    Map = 0;
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
122
    MLPrec = 0;
123

gsell's avatar
gsell committed
124 125
    numBlocks_m = Options::numBlocks;
    recycleBlocks_m = Options::recycleBlocks;
126
    nLHS_m = Options::nLHS;
gsell's avatar
gsell committed
127 128
    SetupMLList();
    SetupBelosList();
129
    problem_ptr = rcp(new Belos::LinearProblem<ST, MV, OP>);
130
    // setup Belos solver
gsell's avatar
gsell committed
131
    if(numBlocks_m == 0 || recycleBlocks_m == 0)
132
        solver_ptr = rcp(new Belos::BlockCGSolMgr<double, MV, OP>());
gsell's avatar
gsell committed
133
    else
134 135
        solver_ptr = rcp(new Belos::RCGSolMgr<double, MV, OP>());
    solver_ptr->setParameters(rcp(&belosList, false));
gsell's avatar
gsell committed
136 137 138 139
    convStatusTest = rcp(new Belos::StatusTestGenResNorm<ST, MV, OP> (tol));
    convStatusTest->defineScaleForm(Belos::NormOfRHS, Belos::TwoNorm);

    //all timers used here
140 141 142 143 144 145 146 147
    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
148 149 150 151 152 153 154 155 156 157
}


MGPoissonSolver::MGPoissonSolver(PartBunch &beam):
    layout_m(&beam.getFieldLayout()), mesh_m(&beam.getMesh()), itsBunch_m(&beam), LHS(0), Comm(Ippl::getComm()) {

    throw OpalException("MGPoissonSolver", "Copy Constructor not implemented yet");
}

MGPoissonSolver::~MGPoissonSolver() {
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
    if (Map) delete Map; Map = 0;
    if (MLPrec) delete MLPrec; MLPrec = 0;
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
    prec_m = Teuchos::null;
    solver_ptr = Teuchos::null;
    problem_ptr = Teuchos::null;
}

void MGPoissonSolver::deletePtr() {
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
    prec_m = Teuchos::null;
gsell's avatar
gsell committed
173 174 175 176 177 178
}

void MGPoissonSolver::computePotential(Field_t &rho, Vector_t hr, double zshift) {
    throw OpalException("MGPoissonSolver", "not implemented yet");
}

179
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
180 181
    if (hasParallelDecompositionChanged_m){
        if (Map !=0 ) delete Map;
gsell's avatar
gsell committed
182
        if(useRCB_m)
183
            redistributeWithRCB(localId);
184
        else if ( bp->getType() == "Geometric" )
185
            IPPLToMap3DGeo(localId);
186
        else if ( bp->getType() == "Elliptic" )
187
            IPPLToMap3D(localId);
188 189

        extrapolateLHS();
190 191
    }
}
gsell's avatar
gsell committed
192

193 194 195 196 197
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.
198
    //we also have to redistribute LHS
199

200 201 202 203 204 205 206 207 208
    if (Teuchos::is_null(LHS)){
        LHS = rcp(new Epetra_Vector(*Map));
        LHS->PutScalar(1.0);
    } else {
        RCP<Epetra_Vector> tmplhs = rcp(new Epetra_Vector(*Map));
        Epetra_Import importer(*Map, LHS->Map());
        tmplhs->Import(*LHS, importer, Add);
        LHS = tmplhs;
    }
gsell's avatar
gsell committed
209

210
    //...and all previously saved LHS
211
    std::deque< Epetra_Vector >::iterator it = OldLHS.begin();
212 213 214 215 216 217 218 219 220 221
    if(OldLHS.size() > 0) {
        int n = OldLHS.size();
        for(int i = 0; i < n; ++i) {
            Epetra_Vector tmplhs = Epetra_Vector(*Map);
            Epetra_Import importer(*Map, it->Map());
            tmplhs.Import(*it, importer, Add);
            *it = tmplhs;
            ++it;
        }
    }
gsell's avatar
gsell committed
222

223 224 225
    // extrapolate last OldLHS.size LHS to get a new start vector
    it = OldLHS.begin();
    if(nLHS_m == 0 || OldLHS.size()==0)
226 227 228
        LHS->PutScalar(1.0);
    else if(OldLHS.size() == 1)
        *LHS = *it;
229 230 231
    else if(OldLHS.size() == 2)
        LHS->Update(2.0, *it++, -1.0, *it, 0.0);
    else if(OldLHS.size() > 2) {
232
        int n = OldLHS.size();
233 234 235
        P = rcp(new Epetra_MultiVector(*Map, nLHS_m, false));
        for(int i = 0; i < n; ++i) {
           *(*P)(i) = *it++;
236
        }
237 238 239 240
        for(int k = 1; k < n; ++k) {
           for(int i = 0; i < n - k; ++i) {
              (*P)(i)->Update(-(i + 1) / (float)k, *(*P)(i + 1), (i + k + 1) / (float)k);
           }
241 242
        }
        *LHS = *(*P)(0);
243 244
     } else
        throw OpalException("MGPoissonSolver", "Invalid number of old LHS: " + OldLHS.size());
gsell's avatar
gsell committed
245 246
}

247

gsell's avatar
gsell committed
248 249 250 251 252
// 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) {
253

gsell's avatar
gsell committed
254 255 256
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
257

258
    //bp->setMinMaxZ(itsBunch_m->get_origin()[2], itsBunch_m->get_maxExtend()[2]);
259 260
    bp->setGlobalMeanR(itsBunch_m->getGlobalMeanR());
    bp->setGlobalToLocalQuaternion(itsBunch_m->getGlobalToLocalQuaternion());
gsell's avatar
gsell committed
261
    bp->setNr(nr_m);
262 263
    NDIndex<3> localId = layout_m->getLocalNDIndex();

264
    IpplTimings::startTimer(FunctionTimer1_m);
265
    if ( bp->getType() == "Geometric" ) {
266
        bp->Compute(hr, localId);
267
    } else if ( bp->getType() == "Elliptic"  )
268
        bp->Compute(hr);
269
    IpplTimings::stopTimer(FunctionTimer1_m);
270 271

    // Define the Map
272
    IpplTimings::startTimer(FunctionTimer2_m);
273
    computeMap(localId);
274
    IpplTimings::stopTimer(FunctionTimer2_m);
gsell's avatar
gsell committed
275

276
    // Allocate the RHS with the new Epetra Map
277
    if (Teuchos::is_null(RHS))
278
      RHS = rcp(new Epetra_Vector(*Map));
279
    RHS->PutScalar(0.0);
280 281 282
    // get charge densities from IPPL field and store in Epetra vector (RHS)
    IpplTimings::startTimer(FunctionTimer3_m);
    int id = 0;
283
    float scaleFactor = itsBunch_m->getdT();
284
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
285
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
286
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
287
                  if (bp->isInside(idx, idy, idz))
288
                    RHS->Values()[id++] = 4*M_PI*rho[idx][idy][idz].get()/scaleFactor;
gsell's avatar
gsell committed
289 290 291
            }
        }
    }
292
    IpplTimings::stopTimer(FunctionTimer3_m);
293

gsell's avatar
gsell committed
294
    // build discretization matrix
295
    IpplTimings::startTimer(FunctionTimer4_m);
296 297
    if (Teuchos::is_null(A))
        A = rcp(new Epetra_CrsMatrix(Copy, *Map,  7, true));
298
    ComputeStencil(hr, RHS);
299
    IpplTimings::stopTimer(FunctionTimer4_m);
gsell's avatar
gsell committed
300 301 302
#ifdef DBG_STENCIL
    EpetraExt::RowMatrixToMatlabFile("DiscrStencil.dat", *A);
#endif
303

304
    IpplTimings::startTimer(FunctionTimer5_m);
305 306 307 308 309 310 311 312
    if ( MLPrec == 0 ){
        MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
    } else if (precmode_m == REUSE_HIERARCHY) {
        MLPrec->ReComputePreconditioner();
    } else if (precmode_m == REUSE_PREC){
    } else if (precmode_m == STD_PREC) {
        delete MLPrec; MLPrec = 0;
        MLPrec = new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m);
313
    }
314
    IpplTimings::stopTimer(FunctionTimer5_m);
315

gsell's avatar
gsell committed
316 317
    // setup preconditioned iterative solver
    // use old LHS solution as initial guess
318
    IpplTimings::startTimer(FunctionTimer6_m);
319 320 321 322 323 324 325 326 327

    problem_ptr->setOperator(A);
    problem_ptr->setLHS(LHS);
    problem_ptr->setRHS(RHS);
    prec_m = Teuchos::rcp ( new Belos::EpetraPrecOp ( rcp(MLPrec,false)));
    problem_ptr->setLeftPrec(prec_m);
    solver_ptr->setProblem( problem_ptr);
    if (!problem_ptr->isProblemSet()) {
        if (problem_ptr->setProblem() == false) {
328
            ERRORMSG("Belos::LinearProblem failed to set up correctly!" << endl);
329
        }
gsell's avatar
gsell committed
330
    }
331
    IpplTimings::stopTimer(FunctionTimer6_m);
gsell's avatar
gsell committed
332 333
    std::ofstream timings;
    char filename[50];
334
    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], Comm.NumProc(), recycleBlocks_m, numBlocks_m, nLHS_m);
gsell's avatar
gsell committed
335 336 337
    if(Comm.MyPID() == 0) timings.open(filename, std::ios::app);
    double time = MPI_Wtime();

338
    IpplTimings::startTimer(FunctionTimer7_m);
339
    	solver_ptr->solve();
340
    IpplTimings::stopTimer(FunctionTimer7_m);
gsell's avatar
gsell committed
341 342 343 344 345 346 347 348

    time = MPI_Wtime() - time;
    double minTime, maxTime, avgTime;
    Comm.MinAll(&time, &minTime, 1);
    Comm.MaxAll(&time, &maxTime, 1);
    Comm.SumAll(&time, &avgTime, 1);
    avgTime /= Comm.NumProc();
    if(Comm.MyPID() == 0) timings <<
349
                                      solver_ptr->getNumIters() << "\t" <<
gsell's avatar
gsell committed
350 351 352 353 354 355
                                      //time <<" "<<
                                      minTime << "\t" <<
                                      maxTime << "\t" <<
                                      avgTime << "\t" <<
                                      numBlocks_m << "\t" <<
                                      recycleBlocks_m << "\t" <<
356
                                      nLHS_m << "\t" <<
gsell's avatar
gsell committed
357 358 359 360 361
                                      //OldLHS.size() <<"\t"<<
                                      std::endl;
    if(Comm.MyPID() == 0) timings.close();

    // Store new LHS in OldLHS
362 363
    if(nLHS_m > 1) OldLHS.push_front(*(LHS.get()));
    if(OldLHS.size() > nLHS_m) OldLHS.pop_back();
gsell's avatar
gsell committed
364 365

    //now transfer solution back to IPPL grid
366 367
    IpplTimings::startTimer(FunctionTimer8_m);
    id = 0;
gsell's avatar
gsell committed
368
    rho = 0.0;
369
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
370
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
371
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
372 373
                  NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
                  if (bp->isInside(idx, idy, idz))
374
                     rho.localElement(l) = LHS->Values()[id++]*scaleFactor;
gsell's avatar
gsell committed
375 376 377
            }
        }
    }
378
    IpplTimings::stopTimer(FunctionTimer8_m);
379 380
    if ( hasParallelDecompositionChanged_m )
        deletePtr();
gsell's avatar
gsell committed
381 382 383
}


384
void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
gsell's avatar
gsell committed
385 386 387

    int numMyGridPoints = 0;

388
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
389
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
390
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
391
                 if (bp->isInside(idx, idy, idz))
gsell's avatar
gsell committed
392
                    numMyGridPoints++;
393
             }
gsell's avatar
gsell committed
394
        }
395
     }
gsell's avatar
gsell committed
396 397 398 399 400 401 402 403 404 405

    Epetra_BlockMap bmap(-1, numMyGridPoints, 1, 0, Comm);
    Teuchos::RCP<const Epetra_MultiVector> coords = Teuchos::rcp(new Epetra_MultiVector(bmap, 3, false));

    double *v;
    int stride, stride2;

    coords->ExtractView(&v, &stride);
    stride2 = 2 * stride;

406
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
407
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
408
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
409 410 411 412
                 if (bp->isInside(idx, idy, idz)) {
                    v[0] = (double)idx;
                    v[stride] = (double)idy;
                    v[stride2] = (double)idz;
gsell's avatar
gsell committed
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441
                    v++;
                }
            }
        }
    }

    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");

    Teuchos::RCP<Isorropia::Epetra::Partitioner> part = Teuchos::rcp(new Isorropia::Epetra::Partitioner(coords, paramlist));
    Isorropia::Epetra::Redistributor rd(part);
    Teuchos::RCP<Epetra_MultiVector> newcoords = rd.redistribute(*coords);

    newcoords->ExtractView(&v, &stride);
    stride2 = 2 * stride;
    numMyGridPoints = 0;
    std::vector<int> MyGlobalElements;
    for(int i = 0; i < newcoords->MyLength(); i++) {
        MyGlobalElements.push_back(bp->getIdx(v[0], v[stride], v[stride2]));
        v++;
        numMyGridPoints++;
    }

    Map = new Epetra_Map(-1, numMyGridPoints, &MyGlobalElements[0], 0, Comm);
}

442
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
gsell's avatar
gsell committed
443
    int NumMyElements = 0;
444
    std::vector<int> MyGlobalElements;
gsell's avatar
gsell committed
445

446
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
447
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
448
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
449
                if (bp->isInside(idx, idy, idz)) {
450
                    MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
451 452 453
                    NumMyElements++;
                }
            }
454 455
	}
    }
gsell's avatar
gsell committed
456 457 458
    Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
}

459 460
void MGPoissonSolver::IPPLToMap3DGeo(NDIndex<3> localId) {
    int NumMyElements = 0;
461
    std::vector<int> MyGlobalElements;
462

463 464 465 466 467 468 469 470 471 472 473
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
                    MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
                    NumMyElements++;
            }
	}
    }
    Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
}

gsell's avatar
gsell committed
474 475
void MGPoissonSolver::ComputeStencil(Vector_t hr, Teuchos::RCP<Epetra_Vector> RHS) {

476 477
    A->PutScalar(0.0);

gsell's avatar
gsell committed
478 479 480
    int NumMyElements = Map->NumMyElements();
    int *MyGlobalElements = Map->MyGlobalElements();

481 482
    std::vector<double> Values(6);
    std::vector<int> Indices(6);
gsell's avatar
gsell committed
483

484
    for(int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
485 486 487

        int NumEntries = 0;

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

491
        bp->getBoundaryStencil(MyGlobalElements[i], WV, EV, SV, NV, FV, BV, CV, scaleFactor);
492
        RHS->Values()[i] *= scaleFactor;
gsell's avatar
gsell committed
493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519

        bp->getNeighbours(MyGlobalElements[i], W, E, S, N, F, B);
        if(E != -1) {
            Indices[NumEntries] = E;
            Values[NumEntries++] = EV;
        }
        if(W != -1) {
            Indices[NumEntries] = W;
            Values[NumEntries++] = WV;
        }
        if(S != -1) {
            Indices[NumEntries] = S;
            Values[NumEntries++] = SV;
        }
        if(N != -1) {
            Indices[NumEntries] = N;
            Values[NumEntries++] = NV;
        }
        if(F != -1) {
            Indices[NumEntries] = F;
            Values[NumEntries++] = FV;
        }
        if(B != -1) {
            Indices[NumEntries] = B;
            Values[NumEntries++] = BV;
        }

520
    	// if matrix has already been filled (FillComplete()) we can only
gsell's avatar
gsell committed
521
        // replace entries
522

523 524 525 526 527
        if(A->Filled()) {
            // off-diagonal entries
            A->ReplaceGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
            // diagonal entry
            A->ReplaceGlobalValues(MyGlobalElements[i], 1, &CV, MyGlobalElements + i);
528
        } else {
gsell's avatar
gsell committed
529 530 531
            // off-diagonal entries
            A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
            // diagonal entry
532
            A->InsertGlobalValues(MyGlobalElements[i], 1, &CV, MyGlobalElements + i);
533
        }
gsell's avatar
gsell committed
534 535 536
    }

    A->FillComplete();
537

gsell's avatar
gsell committed
538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560
    A->OptimizeStorage();
}

void MGPoissonSolver::printLoadBalanceStats() {

    //compute some load balance statistics
    size_t myNumPart = Map->NumMyElements();
    size_t NumPart = Map->NumGlobalElements() * 1.0 / Comm.NumProc();
    double imbalance = 1.0;
    if(myNumPart >= NumPart)
        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());

    avg /= Comm.NumProc();
561 562
    if(Comm.MyPID() == 0) *gmsg << "LBAL min = " << min << ", max = " << max << ", avg = " << avg << endl;
    if(Comm.MyPID() == 0) *gmsg << "min nr gridpoints = " << minn << ", max nr gridpoints = " << maxn << endl;
gsell's avatar
gsell committed
563 564 565 566 567 568
}

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;
569
    return os;
gsell's avatar
gsell committed
570 571
}

572
#endif /* HAVE_SAAMG_SOLVER */