53#include <Tpetra_Import.hpp>
54#include <BelosTpetraAdapter.hpp>
57 #include "TpetraExt_MatrixMatrix.hpp"
60#include "Teuchos_CommandLineProcessor.hpp"
62#include "BelosLinearProblem.hpp"
63#include "BelosRCGSolMgr.hpp"
64#include "BelosBlockCGSolMgr.hpp"
65#include "BelosBiCGStabSolMgr.hpp"
66#include "BelosBlockGmresSolMgr.hpp"
67#include "BelosGCRODRSolMgr.hpp"
69#include <MueLu_CreateTpetraPreconditioner.hpp>
80 std::vector<BoundaryGeometry *> geometries,
86 : isMatrixfilled_m(false)
88 , geometries_m(geometries)
90 , maxiters_m(maxiters)
99 for (
int i = 0; i < 3; i++) {
118 bp_m = std::unique_ptr<IrregularDomain>(
122 bp_m = std::unique_ptr<IrregularDomain>(
131 bp_m = std::unique_ptr<IrregularDomain>(
138 if (localId[0].length() !=
domain_m[0].length() ||
139 localId[1].length() !=
domain_m[1].length()) {
141 "The class ArbitraryDomain only works with parallelization\n"
143 "Please set PARFFTX=FALSE, PARFFTY=FALSE, PARFFTT=TRUE in \n"
144 "the definition of the field solver in the input file.\n");
148 bp_m = std::unique_ptr<IrregularDomain>(
152 map_p = Teuchos::null;
167 if (itsolver ==
"CG") {
177 }
else if (itsolver ==
"BICGSTAB") {
182 }
else if (itsolver ==
"GMRES") {
209 map_p = Teuchos::null;
224 throw OpalException(
"MGPoissonSolver",
"not implemented yet");
243 if (Teuchos::is_null(
LHS)){
248 Tpetra::Import<> importer(
map_p,
LHS->getMap());
249 tmplhs->doImport(*
LHS, importer, Tpetra::CombineMode::ADD);
257 for (
int i = 0; i <
n; ++i) {
259 Tpetra::Import<> importer(
map_p, it->getMap());
260 tmplhs.doImport(*it, importer, Tpetra::CombineMode::ADD);
270 else if (
OldLHS.size() == 1)
272 else if (
OldLHS.size() == 2)
273 LHS->update(2.0, *it++, -1.0, *it, 0.0);
274 else if (
OldLHS.size() > 2) {
277 for (
int i = 0; i <
n; ++i) {
278 *(
P_mp->getVectorNonConst(i)) = *it++;
280 for (
int k = 1; k <
n; ++k) {
281 for (
int i = 0; i <
n - k; ++i) {
282 P_mp->getVectorNonConst(i)->update(-(i + 1) / (float)k, *(
P_mp->getVector(i + 1)), (i + k + 1) / (float)k);
288 "Invalid number of old LHS: " + std::to_string(
OldLHS.size()));
310 bp_m->compute(hr, localId);
321 if (Teuchos::is_null(
RHS))
337 msg <<
"* Node:" <<
Ippl::myNode() <<
", Rho for final element: "
338 << rho[localId[0].last()][localId[1].last()][localId[2].last()].get()
342 msg <<
"* Node:" <<
Ippl::myNode() <<
", Local nx*ny*nz = "
343 << localId[2].last() * localId[0].last() * localId[1].last()
346 <<
", Number of reserved local elements in RHS: "
347 <<
RHS->getLocalLength() <<
endl;
349 <<
", Number of reserved global elements in RHS: "
350 <<
RHS->getGlobalLength() <<
endl;
353 for (
int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
354 for (
int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
355 for (
int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
357 if (
bp_m->isInside(idx, idy, idz))
358 RHS->replaceGlobalValue(
bp_m->getIdx(idx, idy, idz),
368 <<
", Number of Local Inside Points " <<
id <<
endl;
375#if defined(__clang__)
376# pragma clang diagnostic push
377# pragma clang diagnostic warning "-Wdeprecated-declarations"
379 if (Teuchos::is_null(
A))
381#if defined(__clang__)
382# pragma clang diagnostic pop
389 Tpetra::MatrixMarket::Writer<TpetraCrsMatrix_t>::writeSparseFile(
390 "DiscrStencil.dat",
A);
395 if (Teuchos::is_null(
prec_mp)) {
396 Teuchos::RCP<TpetraOperator_t> At = Teuchos::rcp_dynamic_cast<TpetraOperator_t>(
A);
403 MueLu::ReuseTpetraPreconditioner(
A, *
prec_mp);
408 Teuchos::RCP<TpetraOperator_t> At = Teuchos::rcp_dynamic_cast<TpetraOperator_t>(
A);
432 ERRORMSG(
"Belos::LinearProblem failed to set up correctly!" <<
endl);
438 double time = MPI_Wtime();
447 std::ofstream timings;
449 time = MPI_Wtime() - time;
450 double minTime = 0, maxTime = 0, avgTime = 0;
451 reduce(time, minTime, 1, std::less<double>());
452 reduce(time, maxTime, 1, std::greater<double>());
453 reduce(time, avgTime, 1, std::plus<double>());
457 sprintf(filename,
"timing_MX%d_MY%d_MZ%d_nProc%d_recB%d_numB%d_nLHS%d",
461 timings.open(filename, std::ios::app);
462 timings <<
solver_mp->getNumIters() <<
"\t"
485 for (
int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
486 for (
int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
487 for (
int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
489 if (
bp_m->isInside(idx, idy, idz))
507 int NumMyElements = 0;
508 std::vector<TpetraGlobalOrdinal_t> MyGlobalElements;
510 for (
int idz = localId[2].first(); idz <= localId[2].last(); idz++) {
511 for (
int idy = localId[1].first(); idy <= localId[1].last(); idy++) {
512 for (
int idx = localId[0].first(); idx <= localId[0].last(); idx++) {
513 if (
bp_m->isInside(idx, idy, idz)) {
514 MyGlobalElements.push_back(
bp_m->getIdx(idx, idy, idz));
522 map_p = Teuchos::rcp(
new TpetraMap_t(Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid(),
523 &MyGlobalElements[0], NumMyElements, indexbase,
comm_mp));
529 A->setAllToScalar(0.0);
531#if defined(__clang__)
532# pragma clang diagnostic push
533# pragma clang diagnostic warning "-Wdeprecated-declarations"
535 int NumMyElements =
map_p->getNodeNumElements();
536#if defined(__clang__)
537# pragma clang diagnostic pop
540 auto MyGlobalElements =
map_p->getMyGlobalIndices();
542 std::vector<TpetraScalar_t> Values(6);
543 std::vector<TpetraGlobalOrdinal_t> Indices(6);
548 for (
int i = 0 ; i < NumMyElements ; i++) {
552 double scaleFactor=1.0;
554 bp_m->getBoundaryStencil(MyGlobalElements[i], value, scaleFactor);
555 RHS->scale(scaleFactor);
557 bp_m->getNeighbours(MyGlobalElements[i], index);
558 if (index.
east != -1) {
559 Indices[NumEntries] = index.
east;
560 Values[NumEntries++] = value.
east;
562 if (index.
west != -1) {
563 Indices[NumEntries] = index.
west;
564 Values[NumEntries++] = value.
west;
566 if (index.
south != -1) {
567 Indices[NumEntries] = index.
south;
568 Values[NumEntries++] = value.
south;
570 if (index.
north != -1) {
571 Indices[NumEntries] = index.
north;
572 Values[NumEntries++] = value.
north;
574 if (index.
front != -1) {
575 Indices[NumEntries] = index.
front;
576 Values[NumEntries++] = value.
front;
578 if (index.
back != -1) {
579 Indices[NumEntries] = index.
back;
580 Values[NumEntries++] = value.
back;
588 A->replaceGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
590 A->replaceGlobalValues(MyGlobalElements[i], 1, &value.
center, &MyGlobalElements[i]);
593 A->insertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
595 A->insertGlobalValues(MyGlobalElements[i], 1, &value.
center, &MyGlobalElements[i]);
599 RCP<ParameterList_t> params = Teuchos::parameterList();
600 params->set (
"Optimize Storage",
true);
601 A->fillComplete(params);
609#if defined(__clang__)
610# pragma clang diagnostic push
611# pragma clang diagnostic warning "-Wdeprecated-declarations"
613 size_t myNumPart =
map_p->getNodeNumElements();
614 size_t NumPart =
map_p->getGlobalNumElements() * 1.0 /
comm_mp->getSize();
615#if defined(__clang__)
616# pragma clang diagnostic pop
618 double imbalance = 1.0;
619 if (myNumPart >= NumPart)
620 imbalance += (myNumPart - NumPart) / NumPart;
622 imbalance += (NumPart - myNumPart) / NumPart;
624 double max = 0.0,
min = 0.0, avg = 0.0;
625 size_t minn = 0,
maxn = 0;
626 reduce(imbalance,
min, 1, std::less<double>());
627 reduce(imbalance,
max, 1, std::greater<double>());
628 reduce(imbalance, avg, 1, std::plus<double>());
629 reduce(myNumPart, minn, 1, std::less<size_t>());
630 reduce(myNumPart,
maxn, 1, std::greater<size_t>());
633 *
gmsg <<
"LBAL min = " <<
min <<
", max = " <<
max <<
", avg = " << avg <<
endl;
634 *
gmsg <<
"min nr gridpoints = " << minn <<
", max nr gridpoints = " <<
maxn <<
endl;
646 belosList.set(
"Verbosity", Belos::Errors + Belos::Warnings +
647 Belos::TimingDetails + Belos::FinalSummary +
648 Belos::StatusTestDetails);
651 belosList.set(
"Verbosity", Belos::Errors + Belos::Warnings);
664 MueLuList_m.set(
"coarse: max size", coarsest_size);
669 MueLuList_m.set(
"filtered matrix: reuse eigenvalue",
false);
672 MueLuList_m.set(
"repartition: rebalance P and R",
false);
673 MueLuList_m.set(
"repartition: partitioner",
"zoltan2");
674 MueLuList_m.set(
"repartition: min rows per proc", 800);
679 Teuchos::ParameterList smparms;
680 smparms.set(
"chebyshev: degree", 3);
681 smparms.set(
"chebyshev: assume matrix does not change",
false);
682 smparms.set(
"chebyshev: zero starting solution",
true);
683 smparms.set(
"relaxation: sweeps", 3);
713 os <<
"* *************** M G P o i s s o n S o l v e r ************************************ " <<
endl;
714 os <<
"* h " <<
hr_m <<
'\n';
715 os <<
"* ********************************************************************************** " <<
endl;
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Inform & endl(Inform &inf)
Inform & level3(Inform &inf)
constexpr double c
The velocity of light in m/s.
std::string::iterator iterator
int numBlocks
RCG: cycle length.
int nLHS
number of old left hand sides used to extrapolate a new start vector
int recycleBlocks
RCG: number of recycle blocks.
long long getLocalTrackStep() const
virtual Vector_t get_hr() const
int precmode_m
preconditioner mode
Teuchos::RCP< TpetraVector_t > LHS
left hand side of the linear system of equations we solve
void computeMap(NDIndex< 3 > localId)
recomputes the map
Teuchos::MpiComm< int > Comm_t
void setupMueLuList()
Setup the parameters for the SAAMG preconditioner.
IpplTimings::TimerRef FunctionTimer2_m
IpplTimings::TimerRef FunctionTimer4_m
Teuchos::RCP< LinearProblem_t > problem_mp
Tpetra::CrsMatrix TpetraCrsMatrix_t
Inform & print(Inform &os) const
Tpetra::Vector ::scalar_type TpetraScalar_t
void printLoadBalanceStats()
useful load balance information
Tpetra::Vector TpetraVector_t
Vector_t hr_m
mesh spacings in each direction
Teuchos::ParameterList belosList
parameter list for the iterative solver (Belos)
PartBunch * itsBunch_m
PartBunch object.
int recycleBlocks_m
number of vectors in recycle space
int numBlocks_m
maximum number of blocks in Krylov space
IpplTimings::TimerRef FunctionTimer5_m
Tpetra::MultiVector TpetraMultiVector_t
Vektor< int, 3 > nr_m
current number of mesh points in each direction
IpplTimings::TimerRef FunctionTimer3_m
double tol_m
tolerance for the iterative solver
std::deque< TpetraVector_t > OldLHS
void setupBelosList()
Setup the parameters for the Belos iterative solver.
void computePotential(Field_t &rho, Vector_t hr)
IpplTimings::TimerRef FunctionTimer7_m
int maxiters_m
maximal number of iterations for the iterative solver
Teuchos::RCP< MueLuTpetraOperator_t > prec_mp
MueLu preconditioner object.
void IPPLToMap3D(NDIndex< 3 > localId)
void ComputeStencil(Vector_t hr, Teuchos::RCP< TpetraVector_t > RHS)
unsigned int nLHS_m
last N LHS's for extrapolating the new LHS as starting vector
Teuchos::RCP< const Comm_t > comm_mp
communicator used by Trilinos
bool verbose_m
flag specifying if we are verbose
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)
BoundaryGeometry * currentGeometry
holding the currently active geometry
Tpetra::Operator TpetraOperator_t
Vektor< int, 3 > orig_nr_m
global number of mesh points in each direction
Teuchos::RCP< TpetraMap_t > map_p
Map holding the processor distribution of data.
Teuchos::RCP< SolverManager_t > solver_mp
IpplTimings::TimerRef FunctionTimer6_m
Teuchos::RCP< TpetraMultiVector_t > P_mp
IpplTimings::TimerRef FunctionTimer8_m
Teuchos::RCP< TpetraVector_t > RHS
right hand side of our problem
IpplTimings::TimerRef FunctionTimer1_m
Teuchos::ParameterList MueLuList_m
parameter list for the MueLu solver
std::vector< BoundaryGeometry * > geometries_m
container for multiple geometries
Teuchos::RCP< TpetraCrsMatrix_t > A
matrix used in the linear system of equations
std::unique_ptr< IrregularDomain > bp_m
structure that holds boundary points
Topology getTopology() const
std::string getFilename() const
std::vector< unsigned long long > localTimeSteps
Maximal number of timesteps.
static Track * block
The block of track data.
The base class for all OPAL exceptions.
T & localElement(const NDIndex< Dim > &) const
NDIndex< Dim > getLocalNDIndex()
const NDIndex< Dim > & getDomain() const
MFLOAT get_meshSpacing(unsigned d) const
static Communicate * Comm
static TimerRef getTimer(const char *nm)
static void stopTimer(TimerRef t)
static void startTimer(TimerRef t)