MGPoissonSolver.cpp 16.6 KB
Newer Older
1
#ifdef HAVE_SAAMG_SOLVER
gsell's avatar
gsell committed
2
#define DBG_STENCIL
3
#include "Algorithms/PartBunch.h"
gsell's avatar
gsell committed
4 5 6 7 8 9 10 11 12 13 14
#include "MGPoissonSolver.h"
#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"
#include "Utilities/Options.h"
#include <algorithm>
15
#include "omp.h"
gsell's avatar
gsell committed
16 17 18

using Physics::c;

19 20 21 22 23 24 25 26
MGPoissonSolver::MGPoissonSolver(PartBunch &beam,Mesh_t *mesh, FieldLayout_t *fl, std::vector<BoundaryGeometry *> geometries, std::string itsolver, std::string interpl, double tol, int maxiters, std::string precmode):
    itsBunch_m(&beam),
    mesh_m(mesh), 
    layout_m(fl), 
    geometries_m(geometries), 
    tol_m(tol), 
    maxiters_m(maxiters), 
    Comm(Ippl::getComm()) {
gsell's avatar
gsell committed
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
    domain_m = layout_m->getDomain();

    for(int i = 0; i < 3; i++) {
        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;

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

    hasParallelDecompositionChanged_m = true;
    useRCB_m = false;
    if(Ippl::Info->getOutputLevel() > 0)
        verbose_m = true;
    else
        verbose_m = false;

55
    // Find CURRENT geometry
gsell's avatar
gsell committed
56 57
    currentGeometry = geometries_m[0];
    if(currentGeometry->getFilename() == "") {
58
        if(currentGeometry->getTopology() == "ELLIPTIC"){
gsell's avatar
gsell committed
59
            bp = new EllipticDomain(currentGeometry->getA(), currentGeometry->getB(), orig_nr_m, hr_m, interpl);
60
	}
gsell's avatar
gsell committed
61
        else if(currentGeometry->getTopology() == "BOXCORNER") {
62
            bp = new BoxCornerDomain(currentGeometry->getA(), currentGeometry->getB(), currentGeometry->getC(), currentGeometry->getLength(),currentGeometry->getL1(), currentGeometry->getL2(), orig_nr_m, hr_m, interpl);
63
            bp->Compute(itsBunch_m->get_hr());
gsell's avatar
gsell committed
64 65 66 67
        } else {
            ERRORMSG("Geometry not known" << endl);
            exit(1);
        }
68
    } else {
69 70
	localId = layout_m->getLocalNDIndex();
	bp = new ArbitraryDomain(currentGeometry, orig_nr_m, hr_m, interpl);
71
    }
72 73


gsell's avatar
gsell committed
74 75 76 77
    Map = 0;
    A = Teuchos::null;
    LHS = Teuchos::null;
    RHS = Teuchos::null;
78

79
    MLPrec = Teuchos::null;
gsell's avatar
gsell committed
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94
    numBlocks_m = Options::numBlocks;
    recycleBlocks_m = Options::recycleBlocks;
    nLHS = Options::nLHS;
    SetupMLList();
    SetupBelosList();

    // setup Belos solver
    if(numBlocks_m == 0 || recycleBlocks_m == 0)
        solver = rcp(new Belos::BlockCGSolMgr<double, MV, OP>());
    else
        solver = rcp(new Belos::RCGSolMgr<double, MV, OP>());
    convStatusTest = rcp(new Belos::StatusTestGenResNorm<ST, MV, OP> (tol));
    convStatusTest->defineScaleForm(Belos::NormOfRHS, Belos::TwoNorm);

    //all timers used here
95 96 97 98 99 100 101 102
    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
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
}


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() {
    if(Map) delete Map;
}

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

120
void MGPoissonSolver::computeMap(NDIndex<3> localId) {
gsell's avatar
gsell committed
121

122
    if(hasParallelDecompositionChanged_m == true) {
gsell's avatar
gsell committed
123 124 125
        if(Map != 0)
            delete Map;
        if(useRCB_m)
126
            redistributeWithRCB(localId);
gsell's avatar
gsell committed
127
        else
128
            IPPLToMap3D(localId);
129
    }
130
	extrapolateLHS();
131
}
gsell's avatar
gsell committed
132

133 134 135 136 137
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.
gsell's avatar
gsell committed
138

139 140 141 142 143
     if(LHS == Teuchos::null) {
          LHS = rcp(new Epetra_Vector(*Map));
          LHS->PutScalar(1.0);
     } 

144
    std::deque< Epetra_Vector >::iterator it = OldLHS.begin();
gsell's avatar
gsell committed
145

146
    if(nLHS == 0 || OldLHS.size() == 0)
147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
        LHS->PutScalar(1.0);
    else if(OldLHS.size() == 1)
        *LHS = *it;
    else if(OldLHS.size() == 2){
        LHS->Update (2.0, *it++, -1.0, *it, 0.0);
    }
    else if(OldLHS.size() > 0)
    {
        int n = OldLHS.size();
        for(int i=0; i<n; ++i){
            *(*P)(i) = *it++;
        }
        for(int k = 1; k < n; ++k){// x==0 & xi==i+1
            for(int i = 0; i < n-k; ++i){
                (*P)(i)->Update(-(i+1)/(float)k, *(*P)(i+1), (i+k+1)/(float)k);//TODO test
            }
        }
        *LHS = *(*P)(0);
    }
    else
    {
168
        *gmsg << "Invalid number of old LHS: " + OldLHS.size() << endl;
gsell's avatar
gsell committed
169 170 171
    }
}

172

gsell's avatar
gsell committed
173 174 175 176 177 178 179 180
// 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) {
    nr_m[0] = orig_nr_m[0];
    nr_m[1] = orig_nr_m[1];
    nr_m[2] = orig_nr_m[2];
181 182 183
    
    bp->setMinMaxZ(itsBunch_m->get_origin()[2], itsBunch_m->get_maxExtend()[2]);
    
gsell's avatar
gsell committed
184 185
    bp->setNr(nr_m);

186
    localId = layout_m->getLocalNDIndex();
187

188
    IpplTimings::startTimer(FunctionTimer1_m);
189
            bp->Compute(hr, localId, itsBunch_m->getGlobalMeanR(), itsBunch_m->getGlobalToLocalQuaternion());
190 191 192 193 194
    IpplTimings::stopTimer(FunctionTimer1_m);
   
    IpplTimings::startTimer(FunctionTimer2_m);
    	    computeMap(localId);	// Define the Map
    IpplTimings::stopTimer(FunctionTimer2_m);
gsell's avatar
gsell committed
195

196

197 198 199 200 201 202 203 204 205
    // Allocate the RHS with the new Epetra Map
    if(RHS == Teuchos::null) {
    	RHS = rcp(new Epetra_Vector(*Map));
	RHS->PutScalar(0.0);
    }	
 
    // get charge densities from IPPL field and store in Epetra vector (RHS)
    IpplTimings::startTimer(FunctionTimer3_m);
    int id = 0;
gsell's avatar
gsell committed
206

207
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
208
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
209
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
210 211
                 if (bp->isInside(idx, idy, idz)) 
                    RHS->Values()[id++] = rho[idx][idy][idz].get();
gsell's avatar
gsell committed
212 213 214
            }
        }
    }
215
    IpplTimings::stopTimer(FunctionTimer3_m);
216

gsell's avatar
gsell committed
217
    // build discretization matrix
218 219 220
    IpplTimings::startTimer(FunctionTimer4_m);
    if (A == Teuchos::null)
	A = rcp(new Epetra_CrsMatrix(Copy, *Map,  7, true));
221
    ComputeStencil(hr, RHS);
222
    IpplTimings::stopTimer(FunctionTimer4_m);
gsell's avatar
gsell committed
223 224 225 226

#ifdef DBG_STENCIL
    EpetraExt::RowMatrixToMatlabFile("DiscrStencil.dat", *A);
#endif
227

228 229
    IpplTimings::startTimer(FunctionTimer5_m);
    if (MLPrec == Teuchos::null)
230
        MLPrec = rcp(new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m));
231 232 233 234 235 236 237 238 239 240 241 242 243 244 245

    switch(precmode_m) {
        case REUSE_PREC:
            break;

        case REUSE_HIERARCHY:
                MLPrec->ReComputePreconditioner();
            break;

        case STD_PREC:
            MLPrec = rcp(new ML_Epetra::MultiLevelPreconditioner(*A, MLList_m));
            break;
    }   
    IpplTimings::stopTimer(FunctionTimer5_m);
    
gsell's avatar
gsell committed
246 247 248

    // setup preconditioned iterative solver
    // use old LHS solution as initial guess
249
    IpplTimings::startTimer(FunctionTimer6_m);
gsell's avatar
gsell committed
250 251 252 253 254 255
    problem.setOperator(A);
    problem.setLHS(LHS);
    problem.setRHS(RHS);
    prec = rcp(new Belos::EpetraPrecOp(MLPrec));
    problem.setLeftPrec(prec);
    solver->setParameters(rcp(&belosList, false));
256 257 258 259 260
    solver->setProblem(rcp(&problem,false));
    if(!problem.isProblemSet()){
        if (problem.setProblem() == false) {
            std::cout << std::endl << "ERROR:  Belos::LinearProblem failed to set up correctly!" << std::endl;
        }
gsell's avatar
gsell committed
261
    }
262
    IpplTimings::stopTimer(FunctionTimer6_m);
gsell's avatar
gsell committed
263 264 265 266 267 268 269

    std::ofstream timings;
    char filename[50];
    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);
    if(Comm.MyPID() == 0) timings.open(filename, std::ios::app);
    double time = MPI_Wtime();

270
    IpplTimings::startTimer(FunctionTimer7_m);
271
    	solver->solve();
272
    IpplTimings::stopTimer(FunctionTimer7_m);
gsell's avatar
gsell committed
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293

    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 <<
                                      solver->getNumIters() << "\t" <<
                                      //time <<" "<<
                                      minTime << "\t" <<
                                      maxTime << "\t" <<
                                      avgTime << "\t" <<
                                      numBlocks_m << "\t" <<
                                      recycleBlocks_m << "\t" <<
                                      nLHS << "\t" <<
                                      //OldLHS.size() <<"\t"<<
                                      std::endl;
    if(Comm.MyPID() == 0) timings.close();

    // Store new LHS in OldLHS
294
    if(nLHS > 1) OldLHS.push_front(*(LHS.get()));
gsell's avatar
gsell committed
295 296 297
    if(OldLHS.size() > nLHS) OldLHS.pop_back();

    //now transfer solution back to IPPL grid
298 299
    IpplTimings::startTimer(FunctionTimer8_m);
    id = 0;
gsell's avatar
gsell committed
300 301
    rho = 0.0;

302
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
303
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
304
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
305 306 307
                  NDIndex<3> l(Index(idx, idx), Index(idy, idy), Index(idz, idz));
                  if (bp->isInside(idx, idy, idz))
                     rho.localElement(l) = LHS->Values()[id++];
gsell's avatar
gsell committed
308 309 310
            }
        }
    }
311
    IpplTimings::stopTimer(FunctionTimer8_m);
gsell's avatar
gsell committed
312 313 314
}


315
void MGPoissonSolver::redistributeWithRCB(NDIndex<3> localId) {
gsell's avatar
gsell committed
316 317 318

    int numMyGridPoints = 0;

319
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
320
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
321
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
322
                 if (bp->isInside(idx, idy, idz))
gsell's avatar
gsell committed
323
                    numMyGridPoints++;
324
             }
gsell's avatar
gsell committed
325
        }
326
     }
gsell's avatar
gsell committed
327 328 329 330 331 332 333 334 335 336

    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;

337
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
338
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
339
    	    for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
340 341 342 343
                 if (bp->isInside(idx, idy, idz)) {
                    v[0] = (double)idx;
                    v[stride] = (double)idy;
                    v[stride2] = (double)idz;
gsell's avatar
gsell committed
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372
                    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);
}

373
void MGPoissonSolver::IPPLToMap3D(NDIndex<3> localId) {
gsell's avatar
gsell committed
374 375 376 377

    int NumMyElements = 0;
    vector<int> MyGlobalElements;

378
    for (int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
379
        for (int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
380
            for (int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
381 382
                 if (bp->isInside(idx, idy, idz)) {
                    MyGlobalElements.push_back(bp->getIdx(idx, idy, idz));
gsell's avatar
gsell committed
383 384 385
                    NumMyElements++;
                }
            }
386 387
	}
    }
gsell's avatar
gsell committed
388 389 390 391 392
    Map = new Epetra_Map(-1, NumMyElements, &MyGlobalElements[0], 0, Comm);
}

void MGPoissonSolver::ComputeStencil(Vector_t hr, Teuchos::RCP<Epetra_Vector> RHS) {

393 394
    A->PutScalar(0.0);

gsell's avatar
gsell committed
395 396 397 398 399 400
    int NumMyElements = Map->NumMyElements();
    int *MyGlobalElements = Map->MyGlobalElements();

    vector<double> Values(6);
    vector<int> Indices(6);

401
    for(int i = 0 ; i < NumMyElements ; i++) {
gsell's avatar
gsell committed
402 403 404

        int NumEntries = 0;

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

408
        bp->getBoundaryStencil(MyGlobalElements[i], WV, EV, SV, NV, FV, BV, CV, scaleFactor);
409
//        RHS->Values()[i] *= scaleFactor; 
gsell's avatar
gsell committed
410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436

        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;
        }

437
    	// if matrix has already been filled (FillComplete()) we can only
gsell's avatar
gsell committed
438
        // replace entries
439 440 441 442 443 444 445
        
        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);
        } else { 
gsell's avatar
gsell committed
446 447 448
            // off-diagonal entries
            A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
            // diagonal entry
449 450
            A->InsertGlobalValues(MyGlobalElements[i], 1, &CV, MyGlobalElements + i);
        } 
gsell's avatar
gsell committed
451 452 453
    }

    A->FillComplete();
454

gsell's avatar
gsell committed
455
    A->OptimizeStorage();
456

gsell's avatar
gsell committed
457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478
}

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();
479 480
    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
481 482 483 484 485 486
}

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;
487
    return os;	
gsell's avatar
gsell committed
488 489
}

490
#endif /* HAVE_SAAMG_SOLVER */