13 #include <AMReX_ParallelDescriptor.H>
15 #include <boost/filesystem.hpp>
22 const std::string& bsolver,
23 const std::string& prec,
24 const bool& rebalance,
25 const std::string& reuse,
26 const std::string& bcx,
27 const std::string& bcy,
28 const std::string& bcz,
29 const std::string& smoother,
30 const std::size_t& nSweeps,
31 const std::string& interp,
32 const std::string& norm)
34 , comm_mp( new
comm_t( amrex::ParallelDescriptor::Communicator() ) )
46 , fname_m(
OpalData::getInstance()->getInputBasename() + std::string(
".solver"))
47 , flag_m(std::ios::out)
49 node_mp = KokkosClassic::Details::getNode<amr::node_t>();
55 const Boundary bcs[AMREX_SPACEDIM] = {
81 if (boost::filesystem::exists(
fname_m)) {
93 unsigned short baseLevel,
94 unsigned short finestLevel,
106 bool reset = !prevAsGuess;
121 for (
int lev =
nlevel_m - 1; lev > -1; --lev) {
129 for (
int lev = 0; lev <
nlevel_m; ++lev) {
130 int ilev = lbase_m + lev;
132 phi[ilev]->setVal(0.0, phi[ilev]->nGrow());
148 "The number of smoothing sweeps needs to be non-negative!");
156 throw OpalException(
"AmrMultiGrid::setMaxNumberOfIterations()",
157 "The max. number of iterations needs to be positive!");
183 for (
unsigned int i = 0; i < AMREX_SPACEDIM; ++i) {
185 case Boundary::DIRICHLET:
191 case Boundary::PERIODIC:
195 throw OpalException(
"AmrMultiGrid::initPhysicalBoundary_m()",
196 "This type of boundary is not supported");
199 go_t tmp =
bc_m[i]->getNumberOfPoints();
207 const amrex::Vector<AmrGeometry_t>& geom,
220 amrex::Periodicity period(
AmrIntVect_t(D_DECL(0, 0, 0)));
224 for (
int lev = 0; lev <
nlevel_m; ++lev) {
230 rho[ilev]->boxArray(),
231 rho[ilev]->DistributionMap(),
245 mglevel_m[lev]->refmask->setVal(AmrMultiGridLevel_t::Refined::NO, 2);
246 mglevel_m[lev]->refmask->FillBoundary(period);
248 amrex::BoxArray ba =
mglevel_m[lev]->grids;
254 mglevel_m[lev]->crsemask->setVal(AmrMultiGridLevel_t::Refined::NO, 2);
257 for (
int lev = 1; lev <
nlevel_m; ++lev) {
258 mglevel_m[lev]->crsemask->setVal(AmrMultiGridLevel_t::Refined::YES, 0);
261 mglevel_m[lev]->crsemask->setDomainBndry(AmrMultiGridLevel_t::Mask::PHYSBNDRY,
264 mglevel_m[lev]->crsemask->FillBoundary(period);
270 for (
int lev = 0; lev < nlevel_m-1; ++lev) {
272 amrex::BoxArray ba =
mglevel_m[lev]->grids;
274 ba = amrex::intersect(
mglevel_m[lev+1]->grids, ba);
278 amrex::DistributionMapping dmap(ba,
comm_mp->getSize());
280 refined.setVal(AmrMultiGridLevel_t::Refined::YES);
284 mglevel_m[lev]->refmask->copy(refined, 0, 0, 1, 0, 2);
289 mglevel_m[lev]->refmask->setDomainBndry(AmrMultiGridLevel_t::Mask::PHYSBNDRY,
292 mglevel_m[lev]->refmask->FillBoundary(period);
298 for (
int lev = 0; lev <
nlevel_m; ++lev) {
311 for (
int lev = 0; lev <
nlevel_m; ++lev)
319 std::vector<scalar_t> rhsNorms;
320 std::vector<scalar_t> resNorms;
325 [
this](
double& val){ val *=
eps_m; });
344 for (
int lev = 0; lev <
nlevel_m; ++lev) {
359 this->writeResidualNorm_m();
365 return std::accumulate(resNorms.begin(),
367 std::plus<scalar_t>());
372 std::vector<scalar_t>& resNorms)
374 return std::equal(resNorms.begin(), resNorms.end(),
376 std::less<scalar_t>());
381 Teuchos::RCP<vector_t>& r,
382 const Teuchos::RCP<vector_t>& b,
383 const Teuchos::RCP<vector_t>& x)
393 if (
mglevel_m[level]->Bfine_p != Teuchos::null ) {
398 mglevel_m[level]->Awf_p->apply(*x, fine2crse,
403 if (
mglevel_m[level]->Bcrse_p != Teuchos::null ) {
412 Teuchos::RCP<vector_t> uncov_Ax = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p,
true) );
413 mglevel_m[level]->UnCovered_p->apply(fine2crse, *uncov_Ax);
415 Teuchos::RCP<vector_t> uncov_b = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p,
true) );
417 mglevel_m[level]->UnCovered_p->apply(*b, *uncov_b);
422 r->update(1.0, *uncov_b, -1.0, *uncov_Ax, 0.0);
427 Teuchos::RCP<vector_t> Ax = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p,
true) );
430 if (
mglevel_m[level]->Bcrse_p != Teuchos::null ) {
438 r->update(1.0, *b, -1.0, *Ax, 0.0);
450 Teuchos::RCP<vector_t> Ax = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p,
true) );
465 Teuchos::RCP<vector_t> phi_save = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p) );
466 Tpetra::deep_copy(*phi_save, *
mglevel_m[level]->phi_p);
468 mglevel_m[level-1]->error_p->putScalar(0.0);
492 Teuchos::RCP<vector_t> tmp = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p) );
498 Tpetra::deep_copy(*
mglevel_m[level]->residual_p, *tmp);
501 Teuchos::RCP<vector_t> derror = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p,
true) );
507 mglevel_m[level]->error_p->update(1.0, *derror, 1.0);
510 mglevel_m[level]->phi_p->update(1.0, *phi_save, 1.0, *
mglevel_m[level]->error_p, 0.0);
531 Teuchos::RCP<vector_t>& result,
532 const Teuchos::RCP<vector_t>& rhs,
533 const Teuchos::RCP<vector_t>& crs_rhs,
534 const Teuchos::RCP<vector_t>& b)
539 if (
mglevel_m[level]->Bcrse_p != Teuchos::null ) {
540 mglevel_m[level]->Bcrse_p->apply(*crs_rhs, crse2fine);
544 mglevel_m[level]->Anf_p->apply(*rhs, crse2fine,
549 result->update(1.0, *b, -1.0, crse2fine, 0.0);
554 void AmrMultiGrid::writeResidualNorm_m() {
559 out.open(
"residual.dat", std::ios::app);
561 for (
int lev = 0; lev <
nlevel_m; ++lev) {
565 out << std::setw(15) << std::right << tmp;
597 "This type of norm not suppported.");
604 std::vector<scalar_t>& resNorms)
612 if (
comm_mp->getRank() == 0) {
613 out.open(
"residual.dat", std::ios::out);
615 for (
int lev = 0; lev <
nlevel_m; ++lev)
616 out << std::setw(14) << std::right <<
"level" << lev;
621 for (
int lev = 0; lev <
nlevel_m; ++lev) {
631 out << std::setw(15) << std::right << resNorms.back();
648 Teuchos::RCP<vector_t> efield_p = Teuchos::null;
649 for (
int lev = nlevel_m - 1; lev > -1; --lev) {
655 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
656 efield[ilev][d]->setVal(0.0, efield[ilev][d]->nGrow());
657 efield_p->putScalar(0.0);
669 const amrex::Vector<AmrField_u>& phi,
670 const bool& matrices)
698 const amrex::Vector<AmrField_u>& phi,
699 const bool& matrices)
706 D_DECL( invdx[0] * invdx[0],
708 invdx[2] * invdx[2] )
713 mfi.isValid(); ++mfi)
715 const box_t& tbx = mfi.tilebox();
720 const int* lo = tbx.loVect();
721 const int* hi = tbx.hiVect();
723 for (
int i = lo[0]; i <= hi[0]; ++i) {
724 for (
int j = lo[1]; j <= hi[1]; ++j) {
725 #if AMREX_SPACEDIM == 3
726 for (
int k = lo[2]; k <= hi[2]; ++k) {
739 #if AMREX_SPACEDIM == 3
759 const amrex::Vector<AmrField_u>& phi,
760 const bool& matrices)
765 for (
int lev = nlevel_m-1; lev < (int)
smoother_m.size(); ++lev) {
771 for (
int lev = 0; lev <
nlevel_m; ++lev) {
772 this->
open_m(lev, matrices);
782 D_DECL( invdx[0] * invdx[0],
784 invdx[2] * invdx[2] )
788 for (amrex::MFIter mfi(*
mglevel_m[lev]->mask,
true);
789 mfi.isValid(); ++mfi)
791 const box_t& tbx = mfi.tilebox();
798 const int* lo = tbx.loVect();
799 const int* hi = tbx.hiVect();
801 for (
int i = lo[0]; i <= hi[0]; ++i) {
803 for (
int j = lo[1]; j <= hi[1]; ++j) {
805 #if AMREX_SPACEDIM == 3
806 for (
int k = lo[2]; k <= hi[2]; ++k) {
813 D_DECL(ii, jj, kk), rfab);
831 mglevel_m[lev]->rho_p->replaceGlobalValue(gidx, rhofab(iv, 0));
832 mglevel_m[lev]->phi_p->replaceGlobalValue(gidx, pfab(iv, 0));
833 #if AMREX_SPACEDIM == 3
848 if ( matrices && lev >
lbase_m ) {
857 const bool& matrices)
871 Tpetra::StaticProfile) );
880 mglevel_m[level]->Bcrse_p = Teuchos::rcp(
882 Tpetra::StaticProfile) );
894 int nNeighbours = AMREX_D_TERM(2, * 2, * 2);
898 Tpetra::StaticProfile) );
905 nNeighbours = 2 * AMREX_SPACEDIM * AMREX_D_TERM(2, * 2, * 2);
907 mglevel_m[level]->Bfine_p = Teuchos::rcp(
909 Tpetra::StaticProfile) );
917 int nPhysBoundary = 2 * AMREX_SPACEDIM *
nBcPoints_m;
920 int nIntBoundary = AMREX_SPACEDIM *
interface_mp->getNumberOfPoints();
922 int nEntries = (AMREX_SPACEDIM << 1) + 2 + nPhysBoundary + nIntBoundary;
927 Tpetra::StaticProfile) );
934 nEntries = (AMREX_SPACEDIM << 1) + 5 + nPhysBoundary + nIntBoundary;
938 Tpetra::StaticProfile) );
943 mglevel_m[level]->UnCovered_p = Teuchos::rcp(
945 Tpetra::StaticProfile) );
954 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
957 Tpetra::StaticProfile) );
973 const bool& matrices)
997 for (
int d = 0; d < AMREX_SPACEDIM; ++d)
998 mglevel_m[level]->G_p[d]->fillComplete();
1002 mglevel_m[level]->Awf_p->fillComplete();
1004 mglevel_m[level]->UnCovered_p->fillComplete();
1033 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
1034 for (
int shift = -1; shift <= 1; shift += 2) {
1038 switch ( mfab(biv) )
1040 case AmrMultiGridLevel_t::Mask::INTERIOR:
1041 case AmrMultiGridLevel_t::Mask::COVERED:
1043 map[
mglevel_m[level]->serialize(biv)] += invdx2[d];
1046 case AmrMultiGridLevel_t::Mask::BNDRY:
1052 throw OpalException(
"AmrMultiGrid::buildNoFinePoissonMatrix_m()",
1053 "Error in mask for level "
1054 + std::to_string(level) +
"!");
1062 case AmrMultiGridLevel_t::Mask::PHYSBNDRY:
1065 mglevel_m[level]->applyBoundary(biv, d, map,
1070 throw OpalException(
"AmrMultiGrid::buildNoFinePoissonMatrix_m()",
1071 "Error in mask for level "
1072 + std::to_string(level) +
"!");
1078 map[gidx] += AMREX_D_TERM(- 2.0 * invdx2[0],
1084 mglevel_m[level]->Anf_p->insertGlobalValues(gidx,
1104 if ( rfab(iv) == AmrMultiGridLevel_t::Refined::YES )
1127 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
1128 for (
int shift = -1; shift <= 1; shift += 2) {
1132 if ( rfab(biv) != AmrMultiGridLevel_t::Refined::YES )
1137 switch ( mfab(biv) )
1139 case AmrMultiGridLevel_t::Mask::COVERED:
1141 case AmrMultiGridLevel_t::Mask::INTERIOR:
1143 map[
mglevel_m[level]->serialize(biv)] += invdx2[d];
1144 map[gidx] -= invdx2[d];
1147 case AmrMultiGridLevel_t::Mask::BNDRY:
1153 throw OpalException(
"AmrMultiGrid::buildCompositePoissonMatrix_m()",
1154 "Error in mask for level "
1155 + std::to_string(level) +
"!");
1170 map[gidx] -= invdx2[d];
1173 case AmrMultiGridLevel_t::Mask::PHYSBNDRY:
1176 mglevel_m[level]->applyBoundary(biv, d, map,
1180 map[gidx] -= invdx2[d];
1184 throw OpalException(
"AmrMultiGrid::buildCompositePoissonMatrix_m()",
1185 "Error in mask for level "
1186 + std::to_string(level) +
"!");
1207 scalar_t invavg = AMREX_D_PICK(1.0, 0.5, 0.25);
1208 scalar_t value = - invavg * invcdx[d] * invfdx[d];
1210 for (
int d1 = 0; d1 < 2; ++d1) {
1211 #if AMREX_SPACEDIM == 3
1212 for (
int d2 = 0; d2 < 2; ++d2) {
1220 fake[(d+1)%AMREX_SPACEDIM] = d1;
1221 #if AMREX_SPACEDIM == 3
1222 fake[(d+2)%AMREX_SPACEDIM] = d2;
1227 #if AMREX_SPACEDIM == 3
1237 mglevel_m[level]->Awf_p->insertGlobalValues(gidx,
1243 mglevel_m[level]->UnCovered_p->insertGlobalValues(gidx,
1253 D_DECL(
const go_t& ii,
1263 if ( rfab(iv) == AmrMultiGridLevel_t::Refined::NO || level ==
lfine_m )
1275 indices.reserve(2 << (AMREX_SPACEDIM - 1));
1277 values.reserve(2 << (AMREX_SPACEDIM -1));
1280 for (
int iref = 0; iref < 2; ++iref) {
1281 for (
int jref = 0; jref < 2; ++jref) {
1282 #if AMREX_SPACEDIM == 3
1283 for (
int kref = 0; kref < 2; ++kref) {
1285 AmrIntVect_t riv(D_DECL(ii + iref, jj + jref, kk + kref));
1288 values.push_back( AMREX_D_PICK(0.5, 0.25, 0.125) );
1289 #if AMREX_SPACEDIM == 3
1295 mglevel_m[level]->R_p->insertGlobalValues(gidx,
1332 mglevel_m[level]->I_p->insertGlobalValues(gidx,
1362 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
1363 for (
int shift = -1; shift <= 1; shift += 2) {
1368 if ( mfab(niv) == AmrMultiGridLevel_t::Mask::BNDRY )
1376 civ.coarsen(
mglevel_m[level]->refinement());
1382 interface_mp->coarse(civ, map, invdx2[d], d, shift, cfab,
1386 else if ( mfab(niv) == AmrMultiGridLevel_t::Mask::PHYSBNDRY ) {
1387 throw OpalException(
"AmrMultiGrid::buildCrseBoundaryMatrix_m()",
1388 "Fine meshes shouldn't be connected "
1389 "to physical (i.e. mesh) boundary!");
1397 mglevel_m[level]->Bcrse_p->insertGlobalValues(gidx,
1416 if ( rfab(iv) == AmrMultiGridLevel_t::Refined::YES || level ==
lfine_m )
1423 scalar_t invavg = AMREX_D_PICK(1, 0.5, 0.25);
1429 auto fill = [&](
umap_t& map,
1430 D_DECL(
int ii,
int jj,
int kk),
1431 int* begin,
int* end,
int d,
1435 for (
int iref = ii - begin[0]; iref <= ii + end[0]; ++iref) {
1436 for (
int jref = jj - begin[1]; jref <= jj + end[1]; ++jref) {
1437 #if AMREX_SPACEDIM == 3
1438 for (
int kref = kk - begin[2]; kref <= kk + end[2]; ++kref) {
1445 if ( (riv[d] >> 1) == iv[d] ) {
1449 scalar_t value = - invavg * invcdx[d] * invfdx[d];
1454 scalar_t value = invavg * invcdx[d] * invfdx[d];
1455 map[
mglevel_m[level+1]->serialize(riv)] += value;
1457 #if AMREX_SPACEDIM == 3
1471 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
1472 for (
int shift = -1; shift <= 1; shift += 2) {
1475 covered[d] += shift;
1477 if ( rfab(covered) == AmrMultiGridLevel_t::Refined::YES &&
1492 int begin[AMREX_SPACEDIM] = { D_DECL(
int(d == 0),
int(d == 1),
int(d == 2) ) };
1493 int end[AMREX_SPACEDIM] = { D_DECL(
int(d != 0),
int(d != 1),
int(d != 2) ) };
1508 int ii = iv[0] << 1;
1509 int jj = iv[1] << 1;
1510 #if AMREX_SPACEDIM == 3
1511 int kk = iv[2] << 1;
1516 fill(map, D_DECL(ii, jj, kk), &begin[0], &end[0], d, iv, shift, 1.0);
1523 int ii = covered[0] << 1;
1524 int jj = covered[1] << 1;
1525 #if AMREX_SPACEDIM == 3
1526 int kk = covered[2] << 1;
1528 fill(map, D_DECL(ii, jj, kk), &begin[0], &end[0], d, iv, shift, 1.0);
1540 mglevel_m[level]->Bfine_p->insertGlobalValues(gidx,
1578 case AmrMultiGridLevel_t::Mask::INTERIOR:
1580 case AmrMultiGridLevel_t::Mask::COVERED:
1582 map[
mglevel_m[level]->serialize(iv)] -= shift * 0.5 * invdx[dir];
1584 case AmrMultiGridLevel_t::Mask::BNDRY:
1589 throw OpalException(
"AmrMultiGrid::buildGradientMatrix_m()",
1590 "Error in mask for level "
1591 + std::to_string(level) +
"!");
1594 scalar_t value = - shift * 0.5 * invdx[dir];
1600 map[mglevel_m[level]->serialize(tmp)] += 2.0 * value;
1604 map[mglevel_m[level]->serialize(tmp)] -= value;
1607 case AmrMultiGridLevel_t::Mask::PHYSBNDRY:
1611 scalar_t value = - shift * 0.5 * invdx[dir];
1613 mglevel_m[level]->applyBoundary(iv, dir,
1622 for (
int d = 0; d < AMREX_SPACEDIM; ++d) {
1623 for (
int shift = -1; shift <= 1; shift += 2) {
1626 check(niv, mfab, d, shift);
1631 mglevel_m[level]->G_p[d]->insertGlobalValues(gidx,
1642 Teuchos::RCP<vector_t>& mv)
1647 for (amrex::MFIter mfi(mf,
true); mfi.isValid(); ++mfi) {
1648 const amrex::Box& tbx = mfi.tilebox();
1649 const amrex::FArrayBox& fab = mf[mfi];
1651 const int* lo = tbx.loVect();
1652 const int* hi = tbx.hiVect();
1654 for (
int i = lo[0]; i <= hi[0]; ++i) {
1655 for (
int j = lo[1]; j <= hi[1]; ++j) {
1656 #if AMREX_SPACEDIM == 3
1657 for (
int k = lo[2]; k <= hi[2]; ++k) {
1663 mv->replaceGlobalValue(gidx, fab(iv, comp));
1664 #if AMREX_SPACEDIM == 3
1676 const Teuchos::RCP<vector_t>& mv)
1678 Teuchos::ArrayRCP<const amr::scalar_t> data = mv->get1dView();
1680 for (amrex::MFIter mfi(mf,
true); mfi.isValid(); ++mfi) {
1681 const amrex::Box& tbx = mfi.tilebox();
1682 amrex::FArrayBox& fab = mf[mfi];
1684 const int* lo = tbx.loVect();
1685 const int* hi = tbx.hiVect();
1687 for (
int i = lo[0]; i <= hi[0]; ++i) {
1688 for (
int j = lo[1]; j <= hi[1]; ++j) {
1689 #if AMREX_SPACEDIM == 3
1690 for (
int k = lo[2]; k <= hi[2]; ++k) {
1697 fab(iv, comp) = data[lidx];
1699 #if AMREX_SPACEDIM == 3
1714 indices.reserve(map.size());
1715 values.reserve(map.size());
1718 [&](
const std::pair<const go_t, scalar_t>& entry)
1721 if ( entry.first < 0 ) {
1723 "Negative matrix index!");
1727 indices.push_back(entry.first);
1728 values.push_back(entry.second);
1737 Teuchos::RCP<vector_t>&
e,
1738 Teuchos::RCP<vector_t>& r)
1759 Teuchos::RCP<vector_t> tmp = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p) );
1766 mglevel_m[level-1]->residual_p->putScalar(0.0);
1782 fine2crse, Teuchos::NO_TRANS,
1785 if (
mglevel_m[level-1]->Bcrse_p != Teuchos::null ) {
1788 fine2crse, Teuchos::NO_TRANS,
1792 Teuchos::RCP<vector_t> uncoveredRho = Teuchos::rcp(
new vector_t(
mglevel_m[level-1]->map_p,
true) );
1798 Teuchos::RCP<vector_t> tmp2 = Teuchos::rcp(
new vector_t(
mglevel_m[level-1]->map_p,
true) );
1799 mglevel_m[level-1]->UnCovered_p->apply(fine2crse, *tmp2);
1802 mglevel_m[level-1]->residual_p->update(1.0, *uncoveredRho, -1.0, *tmp2, 1.0);
1834 Teuchos::RCP<vector_t> phicrse = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p,
false) );
1839 Teuchos::RCP<vector_t> uncov_phi = Teuchos::rcp(
new vector_t(
mglevel_m[level]->map_p,
true) );
1843 mglevel_m[level]->phi_p->update(1.0, *phicrse, 1.0, *uncov_phi, 0.0);
1852 case Interpolater::TRILINEAR:
1855 case Interpolater::LAGRANGE:
1857 "Not yet implemented.");
1858 case Interpolater::PIECEWISE_CONST:
1863 "No such interpolater available.");
1869 switch ( interface ) {
1870 case Interpolater::TRILINEAR:
1873 case Interpolater::LAGRANGE:
1877 case Interpolater::PIECEWISE_CONST:
1882 "No such interpolater for the coarse-fine interface available.");
1888 const bool& rebalance,
1889 const std::string& reuse)
1893 case BaseSolver::BICGSTAB:
1896 case BaseSolver::MINRES:
1899 case BaseSolver::PCPG:
1902 case BaseSolver::CG:
1905 case BaseSolver::GMRES:
1908 case BaseSolver::STOCHASTIC_CG:
1911 case BaseSolver::RECYCLING_CG:
1914 case BaseSolver::RECYCLING_GMRES:
1918 #ifdef HAVE_AMESOS2_KLU2
1919 case BaseSolver::KLU2:
1923 #if HAVE_AMESOS2_SUPERLU
1924 case BaseSolver::SUPERLU:
1928 #ifdef HAVE_AMESOS2_UMFPACK
1929 case BaseSolver::UMFPACK:
1933 #ifdef HAVE_AMESOS2_PARDISO_MKL
1934 case BaseSolver::PARDISO_MKL:
1938 #ifdef HAVE_AMESOS2_MUMPS
1939 case BaseSolver::MUMPS:
1943 #ifdef HAVE_AMESOS2_LAPACK
1944 case BaseSolver::LAPACK:
1956 "No such bottom solver available.");
1962 const bool& rebalance,
1963 const std::string& reuse)
1986 "No such preconditioner available.");
1993 std::map<std::string, Boundary> map;
1995 map[
"DIRICHLET"] = Boundary::DIRICHLET;
1996 map[
"OPEN"] = Boundary::OPEN;
1997 map[
"PERIODIC"] = Boundary::PERIODIC;
2001 if ( boundary == map.end() )
2002 throw OpalException(
"AmrMultiGrid::convertToEnumBoundary_m()",
2003 "No boundary type '" + bc +
"'.");
2004 return boundary->second;
2009 std::map<std::string, Interpolater> map;
2011 map[
"TRILINEAR"] = Interpolater::TRILINEAR;
2012 map[
"LAGRANGE"] = Interpolater::LAGRANGE;
2013 map[
"PC"] = Interpolater::PIECEWISE_CONST;
2017 if ( interpolater == map.end() )
2018 throw OpalException(
"AmrMultiGrid::convertToEnumInterpolater_m()",
2019 "No interpolater '" + interp +
"'.");
2020 return interpolater->second;
2026 std::map<std::string, BaseSolver> map;
2028 map[
"BICGSTAB"] = BaseSolver::BICGSTAB;
2029 map[
"MINRES"] = BaseSolver::MINRES;
2030 map[
"PCPG"] = BaseSolver::PCPG;
2031 map[
"CG"] = BaseSolver::CG;
2032 map[
"GMRES"] = BaseSolver::GMRES;
2033 map[
"STOCHASTIC_CG"] = BaseSolver::STOCHASTIC_CG;
2034 map[
"RECYCLING_CG"] = BaseSolver::RECYCLING_GMRES;
2035 map[
"RECYCLING_GMRES"] = BaseSolver::RECYCLING_GMRES;
2036 #ifdef HAVE_AMESOS2_KLU2
2037 map[
"KLU2"] = BaseSolver::KLU2;
2039 #ifdef HAVE_AMESOS2_SUPERLU
2040 map[
"SUPERLU"] = BaseSolver::SUPERLU;
2042 #ifdef HAVE_AMESOS2_UMFPACK
2043 map[
"UMFPACK"] = BaseSolver::UMFPACK;
2045 #ifdef HAVE_AMESOS2_PARDISO_MKL
2046 map[
"PARDISO_MKL"] = BaseSolver::PARDISO_MKL;
2048 #ifdef HAVE_AMESOS2_MUMPS
2049 map[
"MUMPS"] = BaseSolver::MUMPS;
2051 #ifdef HAVE_AMESOS2_LAPACK
2052 map[
"LAPACK"] = BaseSolver::LAPACK;
2058 if ( solver == map.end() )
2059 throw OpalException(
"AmrMultiGrid::convertToEnumBaseSolver_m()",
2060 "No bottom solver '" + bsolver +
"'.");
2061 return solver->second;
2067 std::map<std::string, Preconditioner> map;
2077 if ( precond == map.end() )
2078 throw OpalException(
"AmrMultiGrid::convertToEnumPreconditioner_m()",
2079 "No preconditioner '" + prec +
"'.");
2080 return precond->second;
2092 std::map<std::string, Norm> map;
2094 map[
"L1"] = Norm::L1;
2095 map[
"L2"] = Norm::L2;
2096 map[
"LINF"] = Norm::LINF;
2102 if (
n == map.end() )
2104 "No norm '" + norm +
"'.");
2112 std::string dateStr(simtimer.
date());
2113 std::string timeStr(simtimer.
time());
2114 std::string indent(
" ");
2117 outfile <<
"&description\n"
2119 <<
"' " << dateStr <<
"" << timeStr <<
"\",\n"
2120 << indent <<
"contents=\"solver info\"\n"
2122 outfile <<
"¶meter\n"
2123 << indent <<
"name=processors,\n"
2124 << indent <<
"type=long,\n"
2125 << indent <<
"description=\"Number of Cores used\"\n"
2127 outfile <<
"¶meter\n"
2128 << indent <<
"name=revision,\n"
2129 << indent <<
"type=string,\n"
2130 << indent <<
"description=\"git revision of opal\"\n"
2132 outfile <<
"¶meter\n"
2133 << indent <<
"name=flavor,\n"
2134 << indent <<
"type=string,\n"
2135 << indent <<
"description=\"OPAL flavor that wrote file\"\n"
2137 outfile <<
"&column\n"
2138 << indent <<
"name=t,\n"
2139 << indent <<
"type=double,\n"
2140 << indent <<
"units=ns,\n"
2141 << indent <<
"description=\"1 Time\"\n"
2143 outfile <<
"&column\n"
2144 << indent <<
"name=mg_iter,\n"
2145 << indent <<
"type=long,\n"
2146 << indent <<
"units=1,\n"
2147 << indent <<
"description=\"2 Number of Multigrid Iterations\"\n"
2149 outfile <<
"&column\n"
2150 << indent <<
"name=bottom_iter,\n"
2151 << indent <<
"type=long,\n"
2152 << indent <<
"units=1,\n"
2153 << indent <<
"description=\"3 Total Number of Bottom Solver Iterations\"\n"
2155 outfile <<
"&column\n"
2156 << indent <<
"name=regrid,\n"
2157 << indent <<
"type=bool,\n"
2158 << indent <<
"units=1,\n"
2159 << indent <<
"description=\"4 Regrid Step\"\n"
2161 outfile <<
"&column\n"
2162 << indent <<
"name=" +
snorm_m +
",\n"
2163 << indent <<
"type=double,\n"
2164 << indent <<
"units=1,\n"
2165 << indent <<
"description=\"5 Error\"\n"
2168 << indent <<
"mode=ascii,\n"
2169 << indent <<
"no_row_counts=1\n"
2182 unsigned int pwi = 10;
2184 std::ofstream outfile;
2186 if (
comm_mp->getRank() == 0 ) {
2188 outfile.precision(15);
2189 outfile.setf(std::ios::scientific, std::ios::floatfield);
2191 if (
flag_m == std::ios::out ) {
2197 << this->
nIter_m << std::setw(pwi) <<
'\t'
2198 << this->
bIter_m << std::setw(pwi) <<
'\t'
2199 << this->
regrid_m << std::setw(pwi) <<
'\t'
2209 void AmrMultiGrid::initTimer_m() {
2253 os <<
"* ********************* A M R M u l t i G r i d ********************** " <<
endl
2255 <<
"* ******************************************************************** " <<
endl;
The global OPAL structure.
std::vector< std::unique_ptr< AmrMultiGridLevel_t > > mglevel_m
container for levels
static void fillMap(map_t &map)
std::unique_ptr< AmrInterpolater< AmrMultiGridLevel_t > > interface_mp
interpolater for coarse-fine interface
constexpr double e
The value of .
void close_m(const lo_t &level, const bool &matrices)
void buildNoFinePoissonMatrix_m(const lo_t &level, const go_t &gidx, const AmrIntVect_t &iv, const basefab_t &mfab, const scalar_t *invdx2)
MueLuPreconditioner< AmrMultiGridLevel_t > MueLuPreconditioner_t
Norm convertToEnumNorm_m(const std::string &norm)
AmrMultiGrid(AmrBoxLib *itsAmrObject_p, const std::string &bsolver, const std::string &prec, const bool &rebalance, const std::string &reuse, const std::string &bcx, const std::string &bcy, const std::string &bcz, const std::string &smoother, const std::size_t &nSweeps, const std::string &interp, const std::string &norm)
void residual_m(const lo_t &level, Teuchos::RCP< vector_t > &r, const Teuchos::RCP< vector_t > &b, const Teuchos::RCP< vector_t > &x)
void buildDensityVector_m(const lo_t &level, const AmrField_t &rho)
amrex::FArrayBox farraybox_t
Teuchos::RCP< comm_t > comm_mp
communicator
void smooth_m(const lo_t &level, Teuchos::RCP< vector_t > &e, Teuchos::RCP< vector_t > &r)
double getXRangeMin(unsigned short level=0)
The base class for all OPAL exceptions.
void buildRestrictionMatrix_m(const lo_t &level, const go_t &gidx, const AmrIntVect_t &iv, D_DECL(const go_t &ii, const go_t &jj, const go_t &kk), const basefab_t &rfab)
int nBcPoints_m
maximum number of stencils points for BC
scalar_t getLevelResidualNorm(lo_t level)
bool verbose_m
If true, a SDDS file is written.
Smoother
All supported Ifpack2 smoothers.
std::string toUpper(const std::string &str)
void buildPotentialVector_m(const lo_t &level, const AmrField_t &phi)
void initCrseFineInterp_m(const Interpolater &interface)
void serialize(Param_t params, std::ostringstream &os)
serializes params using Boost archive
void setVerbose(bool verbose)
void initGuess_m(bool reset)
Amesos2BottomSolver< AmrMultiGridLevel_t > Amesos2Solver_t
std::ios_base::openmode flag_m
std::ios::out or std::ios::app
void open_m(const lo_t &level, const bool &matrices)
MueLuBottomSolver< AmrMultiGridLevel_t > MueLuSolver_t
BaseSolver convertToEnumBaseSolver_m(const std::string &bsolver)
void initInterpolater_m(const Interpolater &interp)
std::size_t nIter_m
number of iterations till convergence
int nlevel_m
number of levelss
std::size_t nSweeps_m
number of smoothing iterations
AmrMultiGridLevel_t::AmrIntVect_t AmrIntVect_t
void solve(AmrScalarFieldContainer_t &rho, AmrScalarFieldContainer_t &phi, AmrVectorFieldContainer_t &efield, unsigned short baseLevel, unsigned short finestLevel, bool prevAsGuess=true)
std::string date() const
Return date.
double getYRangeMin(unsigned short level=0)
void setNumberOfSweeps(const std::size_t &nSweeps)
static OpalData * getInstance()
#define OPAL_PROJECT_NAME
int lbase_m
base level (currently only 0 supported)
static Smoother convertToEnumSmoother(const std::string &smoother)
void writeSDDSData_m(const scalar_t &error)
double getXRangeMax(unsigned short level=0)
Interpolater convertToEnumInterpolater_m(const std::string &interp)
amrex::FabArray< basefab_t > mask_t
static std::string convertToMueLuReuseOption(const std::string &reuse)
std::string snorm_m
norm for convergence criteria
void buildSingleLevel_m(const amrex::Vector< AmrField_u > &rho, const amrex::Vector< AmrField_u > &phi, const bool &matrices=true)
AmrMultiGridLevel< matrix_t, vector_t > AmrMultiGridLevel_t
void buildCompositePoissonMatrix_m(const lo_t &level, const go_t &gidx, const AmrIntVect_t &iv, const basefab_t &mfab, const basefab_t &rfab, const basefab_t &cfab, const scalar_t *invdx2)
smoothed aggregation multigrid
double getZRangeMax(unsigned short level=0)
std::shared_ptr< bsolver_t > solver_mp
bottom solver
void buildCrseBoundaryMatrix_m(const lo_t &level, const go_t &gidx, const AmrIntVect_t &iv, const basefab_t &mfab, const basefab_t &cfab, const scalar_t *invdx2)
static void startTimer(TimerRef t)
AmrMultiGridLevel_t::coefficients_t coefficients_t
void initPrec_m(const Preconditioner &prec, const bool &rebalance, const std::string &reuse)
std::string getInputFn()
get opals input filename
bool regrid_m
is set to true by itsAmrObject_mp and reset to false by solver
double getYRangeMax(unsigned short level=0)
#define OPAL_PROJECT_VERSION
std::size_t getNumIters()
void map2vector_m(umap_t &map, indices_t &indices, coefficients_t &values)
amr::local_ordinal_t lo_t
BaseSolver
Supported bottom solvers.
std::vector< std::shared_ptr< AmrSmoother > > smoother_m
error smoother
AmrMultiGridLevel_t::AmrField_t AmrField_t
Norm norm_m
norm for convergence criteria (l1, l2, linf)
boundary_t bc_m[AMREX_SPACEDIM]
boundary conditions
void restrict_m(const lo_t &level)
double getZRangeMin(unsigned short level=0)
static std::string convertToMueLuReuseOption(const std::string &reuse)
void initResidual_m(std::vector< scalar_t > &rhsNorms, std::vector< scalar_t > &resNorms)
void averageDown_m(const lo_t &level)
void prolongate_m(const lo_t &level)
Smoother convertToEnumSmoother_m(const std::string &smoother)
std::size_t maxiter_m
maximum number of iterations allowed
void buildInterpolationMatrix_m(const lo_t &level, const go_t &gidx, const AmrIntVect_t &iv, const basefab_t &cfab)
AmrMultiGridLevel_t::umap_t umap_t
Gauss-Seidel block relaxation.
void initPhysicalBoundary_m(const Boundary *bc)
Gauss-Seidel point relaxation.
Interface to Ifpack2 smoothers of the Trilinos package.
std::string time() const
Return time.
Inform & print(Inform &os) const
void initBaseSolver_m(const BaseSolver &solver, const bool &rebalance, const std::string &reuse)
std::string fname_m
SDDS filename.
Norm
Supported convergence criteria.
Ifpack2Preconditioner< AmrMultiGridLevel_t > Ifpack2Preconditioner_t
Boundary convertToEnumBoundary_m(const std::string &bc)
Preconditioner convertToEnumPreconditioner_m(const std::string &prec)
void buildFineBoundaryMatrix_m(const lo_t &level, const go_t &gidx, const AmrIntVect_t &iv, const basefab_t &mfab, const basefab_t &rfab, const basefab_t &cfab)
void setup_m(const amrex::Vector< AmrField_u > &rho, const amrex::Vector< AmrField_u > &phi, const bool &matrices=true)
AmrBoxLib * itsAmrObject_mp
std::string getGitRevision()
Smoother smootherType_m
type of smoother
void residual_no_fine_m(const lo_t &level, Teuchos::RCP< vector_t > &result, const Teuchos::RCP< vector_t > &rhs, const Teuchos::RCP< vector_t > &crs_rhs, const Teuchos::RCP< vector_t > &b)
void writeSDDSHeader_m(std::ofstream &outfile)
Teuchos::RCP< amr::node_t > node_mp
kokkos node
const Vector_t & getMeshScaling() const
Boundary
Supported physical boundaries.
static void fillMap(map_t &map)
void buildGradientMatrix_m(const lo_t &level, const go_t &gidx, const AmrIntVect_t &iv, const basefab_t &mfab, const scalar_t *invdx)
scalar_t evalNorm_m(const Teuchos::RCP< const vector_t > &x)
void buildMultiLevel_m(const amrex::Vector< AmrField_u > &rho, const amrex::Vector< AmrField_u > &phi, const bool &matrices=true)
amr::global_ordinal_t go_t
void trilinos2amrex_m(const lo_t &level, const lo_t &comp, AmrField_t &mf, const Teuchos::RCP< vector_t > &mv)
std::shared_ptr< preconditioner_t > prec_mp
preconditioner for bottom solver
BelosBottomSolver< AmrMultiGridLevel_t > BelosSolver_t
AmrMultiGridLevel_t::indices_t indices_t
static TimerRef getTimer(const char *nm)
void initLevels_m(const amrex::Vector< AmrField_u > &rho, const amrex::Vector< AmrGeometry_t > &geom, bool regrid)
bool for_each(const BareFieldIterator< T, D > &p, SameFieldID s, C)
amr::AmrScalarFieldContainer_t AmrScalarFieldContainer_t
void setMaxNumberOfIterations(const std::size_t &maxiter)
static void stopTimer(TimerRef t)
std::size_t bIter_m
number of iterations of bottom solver
const scalar_t eps_m
rhs scale for convergence
void relax_m(const lo_t &level)
Interpolater
Supported interpolaters for prolongation operation.
void computeEfield_m(AmrVectorFieldContainer_t &efield)
bool isConverged_m(std::vector< scalar_t > &rhsNorms, std::vector< scalar_t > &resNorms)
Inform & endl(Inform &inf)
amrex::BaseFab< int > basefab_t
amr::AmrVectorFieldContainer_t AmrVectorFieldContainer_t
std::unique_ptr< AmrInterpolater< AmrMultiGridLevel_t > > interp_mp
interpolater without coarse-fine interface
void amrex2trilinos_m(const lo_t &level, const lo_t &comp, const AmrField_t &mf, Teuchos::RCP< vector_t > &mv)