3#include <AMReX_Config.H>
16 :
public std::runtime_error
19 using std::runtime_error::runtime_error;
44 template <
typename AMF>
46 RT a_tol_rel,
RT a_tol_abs,
const char* checkpoint_file =
nullptr);
48 template <
typename AMF>
49 RT solve (std::initializer_list<AMF*> a_sol,
50 std::initializer_list<AMF const*> a_rhs,
51 RT a_tol_rel,
RT a_tol_abs,
const char* checkpoint_file =
nullptr);
54 RT a_tol_rel,
RT a_tol_abs);
56 template <
typename AMF>
58 Location a_loc = Location::FaceCenter);
60 template <
typename AMF>
62 Location a_loc = Location::FaceCenter);
67 template <
typename AMF>
69 Location a_loc = Location::FaceCenter);
71 template <
typename AMF>
73 Location a_loc = Location::FaceCenter);
75 template <
typename AMF>
78 Location a_loc = Location::FaceCenter);
80 template <
typename AMF>
82 std::initializer_list<AMF*> a_sol,
83 Location a_loc = Location::FaceCenter);
85 template <
typename AMF>
87 Location a_loc = Location::CellCenter);
89 template <
typename AMF>
90 void getFluxes (std::initializer_list<AMF*> a_flux,
91 Location a_loc = Location::CellCenter);
93 template <
typename AMF>
96 Location a_loc = Location::CellCenter);
98 template <
typename AMF>
99 void getFluxes (std::initializer_list<AMF*> a_flux,
100 std::initializer_list<AMF*> a_sol,
101 Location a_loc = Location::CellCenter);
158#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
169 void setHypreOptionsNamespace(
const std::string& prefix)
noexcept
171 hypre_options_namespace = prefix;
174 void setHypreOldDefault (
bool l)
noexcept {hypre_old_default = l;}
175 void setHypreRelaxType (
int n)
noexcept {hypre_relax_type = n;}
176 void setHypreRelaxOrder (
int n)
noexcept {hypre_relax_order = n;}
177 void setHypreNumSweeps (
int n)
noexcept {hypre_num_sweeps = n;}
178 void setHypreStrongThreshold (Real t)
noexcept {hypre_strong_threshold = t;}
183 template <
typename AMF>
184 void prepareForSolve (Vector<AMF*>
const& a_sol, Vector<AMF const*>
const& a_rhs);
196 void mgVcycle (
int amrlev,
int mglev);
200 void NSolve (MLMGT<MF>& a_solver, MF& a_sol, MF& a_rhs);
220#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
221 template <
class TMF=MF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,
int> = 0>
222 void bottomSolveWithHypre (MF&
x,
const MF&
b);
225#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
226 template <
class TMF=MF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,
int> = 0>
227 void bottomSolveWithPETSc (MF&
x,
const MF&
b);
291#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
296 std::unique_ptr<Hypre> hypre_solver;
297 std::unique_ptr<MLMGBndryT<MF>> hypre_bndry;
298 std::unique_ptr<HypreNodeLap> hypre_node_solver;
300 std::string hypre_options_namespace =
"hypre";
301 bool hypre_old_default =
true;
302 int hypre_relax_type = 6;
303 int hypre_relax_order = 1;
304 int hypre_num_sweeps = 2;
305 Real hypre_strong_threshold = 0.25;
309#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
310 std::unique_ptr<PETScABecLap> petsc_solver;
311 std::unique_ptr<MLMGBndryT<MF>> petsc_bndry;
345 RT a_tol_rel,
RT a_tol_abs,
const char* a_file_name)
const;
349template <
typename MF>
351 : linop(a_lp), ncomp(a_lp.getNComp()), namrlevs(a_lp.NAMRLevels()),
352 finest_amr_lev(a_lp.NAMRLevels()-1)
357template <
typename MF>
358template <
typename AMF>
361 std::initializer_list<AMF const*> a_rhs,
362 RT a_tol_rel,
RT a_tol_abs,
const char* checkpoint_file) ->
RT
366 a_tol_rel, a_tol_abs, checkpoint_file);
369template <
typename MF>
370template <
typename AMF>
373 RT a_tol_rel,
RT a_tol_abs,
const char* checkpoint_file) ->
RT
377 if constexpr (std::is_same<AMF,MultiFab>()) {
378 if (checkpoint_file !=
nullptr) {
379 checkPoint(a_sol, a_rhs, a_tol_rel, a_tol_abs, checkpoint_file);
384 bottom_solver = linop.getDefaultBottomSolver();
387#if (defined(AMREX_USE_HYPRE) || defined(AMREX_USE_PETSC)) && (AMREX_SPACEDIM > 1)
389 int mo = linop.getMaxOrder();
390 if (a_sol[0]->hasEBFabFactory()) {
391 linop.setMaxOrder(2);
393 linop.setMaxOrder(std::min(3,mo));
398 bool is_nsolve = linop.m_parent;
402 RT& composite_norminf = m_final_resnorm0;
405 m_iter_fine_resnorm0.clear();
407 prepareForSolve(a_sol, a_rhs);
409 computeMLResidual(finest_amr_lev);
412 RT resnorm0 = MLResNormInf(finest_amr_lev, local);
413 RT rhsnorm0 = MLRhsNormInf(local);
419 amrex::Print() << print_ident <<
"MLMG: Initial rhs = " << rhsnorm0 <<
"\n"
420 << print_ident <<
"MLMG: Initial residual (resid0) = " << resnorm0 <<
"\n";
424 m_init_resnorm0 = resnorm0;
425 m_rhsnorm0 = rhsnorm0;
428 std::string norm_name;
429 if (always_use_bnorm || rhsnorm0 >= resnorm0) {
433 norm_name =
"resid0";
436 const RT res_target = std::max(a_tol_abs, std::max(a_tol_rel,
RT(1.e-16))*max_norm);
438 if (!is_nsolve && resnorm0 <= res_target) {
439 composite_norminf = resnorm0;
441 amrex::Print() << print_ident <<
"MLMG: No iterations needed\n";
445 bool converged =
false;
447 const int niters = do_fixed_number_of_iters ? do_fixed_number_of_iters : max_iters;
448 for (
int iter = 0; iter < niters; ++iter)
455 computeResidual(finest_amr_lev);
457 if (is_nsolve) {
continue; }
459 RT fine_norminf = ResNormInf(finest_amr_lev);
460 m_iter_fine_resnorm0.push_back(fine_norminf);
461 composite_norminf = fine_norminf;
463 amrex::Print() << print_ident <<
"MLMG: Iteration " << std::setw(3) << iter+1 <<
" Fine resid/"
464 << norm_name <<
" = " << fine_norminf/max_norm <<
"\n";
466 bool fine_converged = (fine_norminf <= res_target);
468 if (namrlevs == 1 && fine_converged) {
470 }
else if (fine_converged) {
472 computeMLResidual(finest_amr_lev-1);
473 RT crse_norminf = MLResNormInf(finest_amr_lev-1);
475 amrex::Print() << print_ident <<
"MLMG: Iteration " << std::setw(3) << iter+1
476 <<
" Crse resid/" << norm_name <<
" = "
477 << crse_norminf/max_norm <<
"\n";
479 converged = (crse_norminf <= res_target);
480 composite_norminf = std::max(fine_norminf, crse_norminf);
487 amrex::Print() << print_ident <<
"MLMG: Final Iter. " << iter+1
488 <<
" resid, resid/" << norm_name <<
" = "
489 << composite_norminf <<
", "
490 << composite_norminf/max_norm <<
"\n";
494 if (composite_norminf >
RT(1.e20)*max_norm)
497 amrex::Print() << print_ident <<
"MLMG: Failing to converge after " << iter+1 <<
" iterations."
498 <<
" resid, resid/" << norm_name <<
" = "
499 << composite_norminf <<
", "
500 << composite_norminf/max_norm <<
"\n";
503 if ( throw_exception ) {
504 throw error(
"MLMG blew up.");
512 if (!converged && do_fixed_number_of_iters == 0) {
514 amrex::Print() << print_ident <<
"MLMG: Failed to converge after " << max_iters <<
" iterations."
515 <<
" resid, resid/" << norm_name <<
" = "
516 << composite_norminf <<
", "
517 << composite_norminf/max_norm <<
"\n";
520 if ( throw_exception ) {
521 throw error(
"MLMG failed to converge.");
532 if (linop.hasHiddenDimension()) {
533 ng_back[linop.hiddenDirection()] = 0;
535 for (
int alev = 0; alev < namrlevs; ++alev)
537 if (!sol_is_alias[alev]) {
538 LocalCopy(*a_sol[alev], sol[alev], 0, 0, ncomp, ng_back);
544 ParallelReduce::Max<double>(timer.data(), timer.size(), 0,
548 amrex::AllPrint() << print_ident <<
"MLMG: Timers: Solve = " << timer[solve_time]
549 <<
" Iter = " << timer[iter_time]
550 <<
" Bottom = " << timer[bottom_time] <<
"\n";
556 return composite_norminf;
559template <
typename MF>
562 RT a_tol_rel,
RT a_tol_abs) ->
RT
565 std::swap(max_precond_iters, do_fixed_number_of_iters);
566 linop.beginPrecondBC();
568 auto r = solve(a_sol, a_rhs, a_tol_rel, a_tol_abs);
570 linop.endPrecondBC();
571 std::swap(max_precond_iters, do_fixed_number_of_iters);
572 precond_mode =
false;
577template <
typename MF>
581 for (
int alev = finest_amr_lev; alev >= 0; --alev) {
582 const MF* crse_bcdata = (alev > 0) ? a_sol[alev-1] :
nullptr;
583 linop.prepareForFluxes(alev, crse_bcdata);
587template <
typename MF>
588template <
typename AMF>
593 for (
int alev = 0; alev <= finest_amr_lev; ++alev) {
594 if constexpr (std::is_same<AMF,MF>()) {
595 linop.compGrad(alev, a_grad_sol[alev], sol[alev], a_loc);
598 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
599 auto const& amf = *(a_grad_sol[alev][idim]);
602 linop.compGrad(alev,
GetArrOfPtrs(grad_sol), sol[alev], a_loc);
603 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
604 LocalCopy(*a_grad_sol[alev][idim], grad_sol[idim], 0, 0, ncomp,
IntVect(0));
610template <
typename MF>
611template <
typename AMF>
618template <
typename MF>
619template <
typename AMF>
624 if (!linop.isCellCentered()) {
625 amrex::Abort(
"Calling wrong getFluxes for nodal solver");
630 if constexpr (std::is_same<AMF,MF>()) {
634 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
635 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
636 auto const& amf = *(a_flux[ilev][idim]);
641 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
642 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
643 LocalCopy(*a_flux[ilev][idim], fluxes[ilev][idim], 0, 0, ncomp,
IntVect(0));
649template <
typename MF>
650template <
typename AMF>
658template <
typename MF>
659template <
typename AMF>
666 if (!linop.isCellCentered()) {
667 amrex::Abort(
"Calling wrong getFluxes for nodal solver");
670 if constexpr (std::is_same<AMF,MF>()) {
671 linop.getFluxes(a_flux, a_sol, a_loc);
674 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
675 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
676 auto const& amf = *(a_flux[ilev][idim]);
682 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
683 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
684 LocalCopy(*a_flux[ilev][idim], fluxes[ilev][idim], 0, 0, ncomp,
IntVect(0));
690template <
typename MF>
691template <
typename AMF>
694 std::initializer_list<AMF*> a_sol,
Location a_loc)
700template <
typename MF>
701template <
typename AMF>
706 if constexpr (std::is_same<AMF,MF>()) {
710 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
711 auto const& amf = *a_flux[ilev];
715 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
721template <
typename MF>
722template <
typename AMF>
729template <
typename MF>
730template <
typename AMF>
737 if constexpr (! std::is_same<AMF,MF>()) {
738 for (
int alev = 0; alev < namrlevs; ++alev) {
743 if (linop.isCellCentered())
746 for (
int alev = 0; alev < namrlevs; ++alev) {
747 for (
int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
750 if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(alev); }
751 ffluxes[alev][idim].define(
amrex::convert(linop.m_grids[alev][mglev],
753 linop.m_dmap[alev][mglev], ncomp, nghost,
MFInfo(),
754 *linop.m_factory[alev][mglev]);
757 if constexpr (std::is_same<AMF,MF>()) {
762 for (
int alev = 0; alev < namrlevs; ++alev) {
771 if constexpr (std::is_same<AMF,MF>()) {
772 linop.getFluxes(a_flux, a_sol);
775 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
776 auto const& amf = *a_flux[ilev];
780 for (
int ilev = 0; ilev < namrlevs; ++ilev) {
787template <
typename MF>
788template <
typename AMF>
791 std::initializer_list<AMF*> a_sol,
Location a_loc)
798template <
typename MF>
802 if (!linop.isCellCentered()) {
810template <
typename MF>
812MLMGT<MF>::getEBFluxes (
const Vector<MF*>& a_eb_flux,
const Vector<MF*>& a_sol)
816 if (!linop.isCellCentered()) {
820 linop.getEBFluxes(a_eb_flux, a_sol);
824template <
typename MF>
832 if (linop.hasHiddenDimension()) { ng_sol[linop.hiddenDirection()] = 0; }
835 sol_is_alias.resize(namrlevs,
true);
836 for (
int alev = 0; alev < namrlevs; ++alev)
838 if (cf_strategy == CFStrategy::ghostnodes ||
nGrowVect(*a_sol[alev]) == ng_sol)
840 sol[alev] = linop.makeAlias(*a_sol[alev]);
841 sol_is_alias[alev] =
true;
845 if (sol_is_alias[alev])
847 sol[alev] = linop.make(alev, 0, ng_sol);
855 const auto& amrrr = linop.AMRRefRatio();
857 for (
int alev = finest_amr_lev; alev >= 0; --alev) {
858 const MF* crse_bcdata = (alev > 0) ? &(sol[alev-1]) :
nullptr;
859 const MF* prhs = a_rhs[alev];
860#if (AMREX_SPACEDIM != 3)
861 int nghost = (cf_strategy == CFStrategy::ghostnodes) ? linop.getNGrow(alev) : 0;
863 MFInfo(), *linop.Factory(alev));
865 linop.applyMetricTerm(alev, 0, rhstmp);
866 linop.unimposeNeumannBC(alev, rhstmp);
867 linop.applyInhomogNeumannTerm(alev, rhstmp);
870 linop.solutionResidual(alev, *a_res[alev], sol[alev], *prhs, crse_bcdata);
871 if (alev < finest_amr_lev) {
872 linop.reflux(alev, *a_res[alev], sol[alev], *prhs,
873 *a_res[alev+1], sol[alev+1], *a_rhs[alev+1]);
874 if (linop.isCellCentered()) {
878 average_down(*a_res[alev+1], *a_res[alev], 0, ncomp, amrrr[alev]);
885#if (AMREX_SPACEDIM != 3)
886 for (
int alev = 0; alev <= finest_amr_lev; ++alev) {
887 linop.unapplyMetricTerm(alev, 0, *a_res[alev]);
892template <
typename MF>
903 if (linop.hasHiddenDimension()) { ng_sol[linop.hiddenDirection()] = 0; }
905 for (
int alev = 0; alev < namrlevs; ++alev)
907 if (cf_strategy == CFStrategy::ghostnodes)
909 nghost = linop.getNGrow(alev);
910 in[alev] = a_in[alev];
912 else if (
nGrowVect(*a_in[alev]) == ng_sol)
914 in[alev] = a_in[alev];
919 if (cf_strategy == CFStrategy::ghostnodes) { ng =
IntVect(nghost); }
920 in_raii[alev] = linop.make(alev, 0, ng);
922 in[alev] = &(in_raii[alev]);
924 rh[alev] = linop.make(alev, 0,
IntVect(nghost));
930 for (
int alev = 0; alev < namrlevs; ++alev) {
931 linop.applyInhomogNeumannTerm(alev, rh[alev]);
934 const auto& amrrr = linop.AMRRefRatio();
936 for (
int alev = finest_amr_lev; alev >= 0; --alev) {
937 const MF* crse_bcdata = (alev > 0) ? in[alev-1] :
nullptr;
938 linop.solutionResidual(alev, *out[alev], *in[alev], rh[alev], crse_bcdata);
939 if (alev < finest_amr_lev) {
940 linop.reflux(alev, *out[alev], *in[alev], rh[alev],
941 *out[alev+1], *in[alev+1], rh[alev+1]);
942 if (linop.isCellCentered()) {
943 if constexpr (IsMultiFabLike_v<MF>) {
950 amrex::Abort(
"MLMG: TODO average_down for non-MultiFab");
956#if (AMREX_SPACEDIM != 3)
957 for (
int alev = 0; alev <= finest_amr_lev; ++alev) {
958 linop.unapplyMetricTerm(alev, 0, *out[alev]);
962 for (
int alev = 0; alev <= finest_amr_lev; ++alev) {
963 if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(alev); }
964 Scale(*out[alev],
RT(-1), 0,
nComp(*out[alev]), nghost);
968template <
typename MF>
973 linop.beginPrecondBC();
975 linop.endPrecondBC();
976 precond_mode =
false;
979template <
typename MF>
980template <
typename AMF>
989 timer.assign(ntimers, 0.0);
993 if (linop.hasHiddenDimension()) { ng_sol[linop.hiddenDirection()] = 0; }
995 if (!linop_prepared) {
996 linop.prepareForSolve();
997 linop_prepared =
true;
998 }
else if (linop.needsUpdate()) {
1001#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
1002 hypre_solver.reset();
1003 hypre_bndry.reset();
1004 hypre_node_solver.reset();
1007#ifdef AMREX_USE_PETSC
1008 petsc_solver.reset();
1009 petsc_bndry.reset();
1014 sol_is_alias.resize(namrlevs,
false);
1015 for (
int alev = 0; alev < namrlevs; ++alev)
1017 if (cf_strategy == CFStrategy::ghostnodes)
1019 if constexpr (std::is_same<AMF,MF>()) {
1020 sol[alev] = linop.makeAlias(*a_sol[alev]);
1021 sol_is_alias[alev] =
true;
1023 amrex::Abort(
"Type conversion not supported for CFStrategy::ghostnodes");
1028 if (
nGrowVect(*a_sol[alev]) == ng_sol) {
1029 if constexpr (std::is_same<AMF,MF>()) {
1030 sol[alev] = linop.makeAlias(*a_sol[alev]);
1031 sol_is_alias[alev] =
true;
1034 if (!sol_is_alias[alev]) {
1035 if (!solve_called) {
1036 sol[alev] = linop.make(alev, 0, ng_sol);
1044 rhs.resize(namrlevs);
1045 for (
int alev = 0; alev < namrlevs; ++alev)
1047 if (cf_strategy == CFStrategy::ghostnodes) { ng_rhs =
IntVect(linop.getNGrow(alev)); }
1048 if (!solve_called) {
1049 rhs[alev] = linop.make(alev, 0, ng_rhs);
1051 LocalCopy(rhs[alev], *a_rhs[alev], 0, 0, ncomp, ng_rhs);
1052 linop.applyMetricTerm(alev, 0, rhs[alev]);
1053 linop.unimposeNeumannBC(alev, rhs[alev]);
1054 linop.applyInhomogNeumannTerm(alev, rhs[alev]);
1055 linop.applyOverset(alev, rhs[alev]);
1056 if ( ! precond_mode) {
1057 bool r = linop.scaleRHS(alev, &(rhs[alev]));
1063 if (factory && !factory->isAllRegular()) {
1064 if constexpr (std::is_same<MF,MultiFab>()) {
1068 amrex::Abort(
"TODO: MLMG with EB only works with MultiFab");
1074 for (
int falev = finest_amr_lev; falev > 0; --falev)
1076 linop.averageDownSolutionRHS(falev-1, sol[falev-1], rhs[falev-1], sol[falev], rhs[falev]);
1080 if (linop.isSingular(0) && linop.getEnforceSingularSolvable())
1085 IntVect ng = linop.getNGrowVectRestriction();
1086 if (cf_strategy == CFStrategy::ghostnodes) { ng = ng_rhs; }
1087 if (!solve_called) {
1088 linop.make(res, ng);
1089 linop.make(rescor, ng);
1091 for (
int alev = 0; alev <= finest_amr_lev; ++alev)
1093 const int nmglevs = linop.NMGLevels(alev);
1094 for (
int mglev = 0; mglev < nmglevs; ++mglev)
1096 setVal(res [alev][mglev],
RT(0.0));
1097 setVal(rescor[alev][mglev],
RT(0.0));
1101 if (cf_strategy != CFStrategy::ghostnodes) { ng = ng_sol; }
1103 for (
int alev = 0; alev <= finest_amr_lev; ++alev)
1105 const int nmglevs = linop.NMGLevels(alev);
1106 cor[alev].resize(nmglevs);
1107 for (
int mglev = 0; mglev < nmglevs; ++mglev)
1109 if (!solve_called) {
1111 if (cf_strategy == CFStrategy::ghostnodes) { _ng=
IntVect(linop.getNGrow(alev,mglev)); }
1112 cor[alev][mglev] = linop.make(alev, mglev, _ng);
1118 cor_hold.resize(std::max(namrlevs-1,1));
1121 const int nmglevs = linop.NMGLevels(alev);
1122 cor_hold[alev].resize(nmglevs);
1123 for (
int mglev = 0; mglev < nmglevs-1; ++mglev)
1125 if (!solve_called) {
1127 if (cf_strategy == CFStrategy::ghostnodes) { _ng=
IntVect(linop.getNGrow(alev,mglev)); }
1128 cor_hold[alev][mglev] = linop.make(alev, mglev, _ng);
1130 setVal(cor_hold[alev][mglev],
RT(0.0));
1133 for (
int alev = 1; alev < finest_amr_lev; ++alev)
1135 cor_hold[alev].resize(1);
1136 if (!solve_called) {
1138 if (cf_strategy == CFStrategy::ghostnodes) { _ng=
IntVect(linop.getNGrow(alev)); }
1139 cor_hold[alev][0] = linop.make(alev, 0, _ng);
1141 setVal(cor_hold[alev][0],
RT(0.0));
1145 || !linop.supportNSolve())
1150 if (do_nsolve && ns_linop ==
nullptr)
1156 amrex::Print() << print_ident <<
"MLMG: # of AMR levels: " << namrlevs <<
"\n"
1157 << print_ident <<
" # of MG levels on the coarsest AMR level: " << linop.NMGLevels(0)
1160 amrex::Print() << print_ident <<
" # of MG levels in N-Solve: " << ns_linop->NMGLevels(0) <<
"\n"
1161 << print_ident <<
" # of grids in N-Solve: " << ns_linop->m_grids[0][0].size() <<
"\n";
1166template <
typename MF>
1170 if (!linop_prepared) {
1171 linop.prepareForSolve();
1172 linop_prepared =
true;
1173 }
else if (linop.needsUpdate()) {
1178template <
typename MF>
1183 linop.preparePrecond();
1186template <
typename MF>
1190 if constexpr (IsMultiFabLike_v<MF>) {
1191 ns_linop = linop.makeNLinOp(nsolve_grid_size);
1194 if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(); }
1196 const BoxArray& ba = (*ns_linop).m_grids[0][0];
1200 if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; }
1201 ns_sol = std::make_unique<MF>(ba, dm, ncomp, ng,
MFInfo(), *(ns_linop->Factory(0,0)));
1203 if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; }
1204 ns_rhs = std::make_unique<MF>(ba, dm, ncomp, ng,
MFInfo(), *(ns_linop->Factory(0,0)));
1208 ns_linop->setLevelBC(0, ns_sol.get());
1210 ns_mlmg = std::make_unique<MLMGT<MF>>(*ns_linop);
1211 ns_mlmg->setVerbose(0);
1212 ns_mlmg->setFixedIter(1);
1213 ns_mlmg->setMaxFmgIter(20);
1220template <
typename MF>
1225 for (
int alev = finest_amr_lev; alev > 0; --alev)
1230 if (cf_strategy == CFStrategy::ghostnodes) { nghost =
IntVect(linop.getNGrow(alev)); }
1231 LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost);
1234 computeResWithCrseSolFineCor(alev-1,alev);
1236 if (alev != finest_amr_lev) {
1237 std::swap(cor_hold[alev][0], cor[alev][0]);
1244 if (linop.isSingular(0) && linop.getEnforceSingularSolvable())
1246 makeSolvable(0,0,res[0][0]);
1249 if (iter < max_fmg_iters) {
1256 if (cf_strategy == CFStrategy::ghostnodes) { nghost =
IntVect(linop.getNGrow(0)); }
1257 LocalAdd(sol[0], cor[0][0], 0, 0, ncomp, nghost);
1260 for (
int alev = 1; alev <= finest_amr_lev; ++alev)
1263 interpCorrection(alev);
1266 if (cf_strategy == CFStrategy::ghostnodes) { nghost =
IntVect(linop.getNGrow(alev)); }
1267 LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost);
1269 if (alev != finest_amr_lev) {
1270 LocalAdd(cor_hold[alev][0], cor[alev][0], 0, 0, ncomp, nghost);
1274 computeResWithCrseCorFineCor(alev);
1278 LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost);
1280 if (alev != finest_amr_lev) {
1281 LocalAdd(cor[alev][0], cor_hold[alev][0], 0, 0, ncomp, nghost);
1285 linop.averageDownAndSync(sol);
1288template <
typename MF>
1293 const int mglev = 0;
1294 mgVcycle(amrlev, mglev);
1299template <
typename MF>
1305 const int mglev_bottom = linop.NMGLevels(amrlev) - 1;
1307 for (
int mglev = mglev_top; mglev < mglev_bottom; ++mglev)
1309 BL_PROFILE_VAR(
"MLMG::mgVcycle_down::"+std::to_string(mglev), blp_mgv_down_lev);
1314 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev
1315 <<
" DN: Norm before smooth " <<
norm <<
"\n";
1318 setVal(cor[amrlev][mglev],
RT(0.0));
1319 bool skip_fillboundary =
true;
1320 for (
int i = 0; i < nu1; ++i) {
1321 linop.smooth(amrlev, mglev, cor[amrlev][mglev], res[amrlev][mglev], skip_fillboundary);
1322 skip_fillboundary =
false;
1326 computeResOfCorrection(amrlev, mglev);
1331 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev
1332 <<
" DN: Norm after smooth " <<
norm <<
"\n";
1336 linop.restriction(amrlev, mglev+1, res[amrlev][mglev+1], rescor[amrlev][mglev]);
1345 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev_bottom
1346 <<
" DN: Norm before bottom " <<
norm <<
"\n";
1351 computeResOfCorrection(amrlev, mglev_bottom);
1353 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev_bottom
1354 <<
" UP: Norm after bottom " <<
norm <<
"\n";
1362 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev_bottom
1363 <<
" Norm before smooth " <<
norm <<
"\n";
1365 setVal(cor[amrlev][mglev_bottom],
RT(0.0));
1366 bool skip_fillboundary =
true;
1367 for (
int i = 0; i < nu1; ++i) {
1368 linop.smooth(amrlev, mglev_bottom, cor[amrlev][mglev_bottom],
1369 res[amrlev][mglev_bottom], skip_fillboundary);
1370 skip_fillboundary =
false;
1374 computeResOfCorrection(amrlev, mglev_bottom);
1376 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev_bottom
1377 <<
" Norm after smooth " <<
norm <<
"\n";
1382 for (
int mglev = mglev_bottom-1; mglev >= mglev_top; --mglev)
1384 BL_PROFILE_VAR(
"MLMG::mgVcycle_up::"+std::to_string(mglev), blp_mgv_up_lev);
1386 addInterpCorrection(amrlev, mglev);
1389 computeResOfCorrection(amrlev, mglev);
1391 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev
1392 <<
" UP: Norm before smooth " <<
norm <<
"\n";
1394 for (
int i = 0; i < nu2; ++i) {
1395 linop.smooth(amrlev, mglev, cor[amrlev][mglev], res[amrlev][mglev]);
1398 if (cf_strategy == CFStrategy::ghostnodes) { computeResOfCorrection(amrlev, mglev); }
1402 computeResOfCorrection(amrlev, mglev);
1404 amrex::Print() << print_ident <<
"AT LEVEL " << amrlev <<
" " << mglev
1405 <<
" UP: Norm after smooth " <<
norm <<
"\n";
1413template <
typename MF>
1423 const int amrlev = 0;
1424 const int mg_bottom_lev = linop.NMGLevels(amrlev) - 1;
1426 if (cf_strategy == CFStrategy::ghostnodes) { nghost =
IntVect(linop.getNGrow(amrlev)); }
1428 for (
int mglev = 1; mglev <= mg_bottom_lev; ++mglev)
1430 linop.avgDownResMG(mglev, res[amrlev][mglev], res[amrlev][mglev-1]);
1435 for (
int mglev = mg_bottom_lev-1; mglev >= 0; --mglev)
1438 interpCorrection(amrlev, mglev);
1441 computeResOfCorrection(amrlev, mglev);
1443 LocalCopy(res[amrlev][mglev], rescor[amrlev][mglev], 0, 0, ncomp, nghost);
1446 std::swap(cor[amrlev][mglev], cor_hold[amrlev][mglev]);
1447 mgVcycle(amrlev, mglev);
1448 LocalAdd(cor[amrlev][mglev], cor_hold[amrlev][mglev], 0, 0, ncomp, nghost);
1455template <
typename MF>
1461 NSolve(*ns_mlmg, *ns_sol, *ns_rhs);
1465 actualBottomSolve();
1469template <
typename MF>
1477 MF
const& res_bottom = res[0].back();
1488 RT(-1.0),
RT(-1.0));
1490 linop.copyNSolveSolution(cor[0].back(), a_sol);
1493template <
typename MF>
1499 if (!linop.isBottomActive()) {
return; }
1505 const int amrlev = 0;
1506 const int mglev = linop.NMGLevels(amrlev) - 1;
1507 auto&
x = cor[amrlev][mglev];
1508 auto&
b = res[amrlev][mglev];
1514 bool skip_fillboundary =
true;
1515 for (
int i = 0; i < nuf; ++i) {
1516 linop.smooth(amrlev, mglev,
x,
b, skip_fillboundary);
1517 skip_fillboundary =
false;
1524 if (linop.isBottomSingular() && linop.getEnforceSingularSolvable())
1527 raii_b = linop.make(amrlev, mglev, ng);
1531 makeSolvable(amrlev,mglev,*bottom_b);
1536#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
1537 if constexpr (std::is_same<MF,MultiFab>()) {
1538 bottomSolveWithHypre(
x, *bottom_b);
1542 amrex::Abort(
"Using Hypre as bottom solver not supported in this case");
1547#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
1548 if constexpr (std::is_same<MF,MultiFab>()) {
1549 bottomSolveWithPETSc(
x, *bottom_b);
1553 amrex::Abort(
"Using PETSc as bottom solver not supported in this case");
1566 int ret = bottomSolveWithCG(
x, *bottom_b, cg_type);
1576 setVal(cor[amrlev][mglev],
RT(0.0));
1577 ret = bottomSolveWithCG(
x, *bottom_b, cg_type);
1588 if (ret != 0 && ret != 9) {
1589 setVal(cor[amrlev][mglev],
RT(0.0));
1591 const int n = (ret==0) ? nub : nuf;
1592 for (
int i = 0; i < n; ++i) {
1593 linop.smooth(amrlev, mglev,
x,
b);
1600 if (! timer.empty()) {
1605template <
typename MF>
1615 if (cf_strategy == CFStrategy::ghostnodes) { cg_solver.
setNGhost(linop.getNGrow()); }
1617 int ret = cg_solver.
solve(
x,
b, bottom_reltol, bottom_abstol);
1618 if (ret != 0 &&
verbose > 1) {
1619 amrex::Print() << print_ident <<
"MLMG: Bottom solve failed.\n";
1626template <
typename MF>
1632 const int mglev = 0;
1633 for (
int alev = amrlevmax; alev >= 0; --alev) {
1634 const MF* crse_bcdata = (alev > 0) ? &(sol[alev-1]) :
nullptr;
1635 linop.solutionResidual(alev, res[alev][mglev], sol[alev], rhs[alev], crse_bcdata);
1636 if (alev < finest_amr_lev) {
1637 linop.reflux(alev, res[alev][mglev], sol[alev], rhs[alev],
1638 res[alev+1][mglev], sol[alev+1], rhs[alev+1]);
1644template <
typename MF>
1649 const MF* crse_bcdata = (alev > 0) ? &(sol[alev-1]) :
nullptr;
1650 linop.solutionResidual(alev, res[alev][0], sol[alev], rhs[alev], crse_bcdata);
1654template <
typename MF>
1658 BL_PROFILE(
"MLMG::computeResWithCrseSolFineCor()");
1661 if (cf_strategy == CFStrategy::ghostnodes) {
1662 nghost =
IntVect(std::min(linop.getNGrow(falev),linop.getNGrow(calev)));
1665 MF& crse_sol = sol[calev];
1666 const MF& crse_rhs = rhs[calev];
1667 MF& crse_res = res[calev][0];
1669 MF& fine_sol = sol[falev];
1670 const MF& fine_rhs = rhs[falev];
1671 MF& fine_cor = cor[falev][0];
1672 MF& fine_res = res[falev][0];
1673 MF& fine_rescor = rescor[falev][0];
1675 const MF* crse_bcdata = (calev > 0) ? &(sol[calev-1]) :
nullptr;
1676 linop.solutionResidual(calev, crse_res, crse_sol, crse_rhs, crse_bcdata);
1678 linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res, BCMode::Homogeneous);
1679 LocalCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost);
1681 linop.reflux(calev, crse_res, crse_sol, crse_rhs, fine_res, fine_sol, fine_rhs);
1683 linop.avgDownResAmr(calev, crse_res, fine_res);
1687template <
typename MF>
1691 BL_PROFILE(
"MLMG::computeResWithCrseCorFineCor()");
1694 if (cf_strategy == CFStrategy::ghostnodes) {
1695 nghost =
IntVect(linop.getNGrow(falev));
1698 const MF& crse_cor = cor[falev-1][0];
1700 MF& fine_cor = cor [falev][0];
1701 MF& fine_res = res [falev][0];
1702 MF& fine_rescor = rescor[falev][0];
1705 linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res,
1706 BCMode::Inhomogeneous, &crse_cor);
1707 LocalCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost);
1711template <
typename MF>
1718 if (cf_strategy == CFStrategy::ghostnodes) {
1719 nghost =
IntVect(linop.getNGrow(alev));
1722 MF
const& crse_cor = cor[alev-1][0];
1723 MF & fine_cor = cor[alev ][0];
1725 const Geometry& crse_geom = linop.Geom(alev-1,0);
1728 int ng_dst = linop.isCellCentered() ? 1 : 0;
1729 if (cf_strategy == CFStrategy::ghostnodes)
1731 ng_src = linop.getNGrow(alev-1);
1732 ng_dst = linop.getNGrow(alev-1);
1735 MF cfine = linop.makeCoarseAmr(alev,
IntVect(ng_dst));
1740 linop.interpolationAmr(alev, fine_cor, cfine, nghost);
1746template <
typename MF>
1752 MF& crse_cor = cor[alev][mglev+1];
1753 MF& fine_cor = cor[alev][mglev ];
1754 linop.interpAssign(alev, mglev, fine_cor, crse_cor);
1758template <
typename MF>
1764 const MF& crse_cor = cor[alev][mglev+1];
1765 MF& fine_cor = cor[alev][mglev ];
1770 if (linop.isMFIterSafe(alev, mglev, mglev+1))
1776 cfine = linop.makeCoarseMG(alev, mglev,
IntVect(0));
1781 linop.interpolation(alev, mglev, fine_cor, *cmf);
1788template <
typename MF>
1793 MF &
x = cor[amrlev][mglev];
1794 const MF&
b = res[amrlev][mglev];
1795 MF &
r = rescor[amrlev][mglev];
1796 linop.correctionResidual(amrlev, mglev,
r,
x,
b, BCMode::Homogeneous);
1800template <
typename MF>
1805 return linop.normInf(alev, res[alev][0], local);
1809template <
typename MF>
1815 for (
int alev = 0; alev <= alevmax; ++alev)
1817 r = std::max(
r, ResNormInf(alev,
true));
1824template <
typename MF>
1830 for (
int alev = 0; alev <= finest_amr_lev; ++alev) {
1831 auto t = linop.normInf(alev, rhs[alev],
true);
1838template <
typename MF>
1842 auto const&
offset = linop.getSolvabilityOffset(0, 0, rhs[0]);
1844 for (
int c = 0; c < ncomp; ++c) {
1845 amrex::Print() << print_ident <<
"MLMG: Subtracting " <<
offset[c] <<
" from rhs component "
1849 for (
int alev = 0; alev < namrlevs; ++alev) {
1850 linop.fixSolvabilityByOffset(alev, 0, rhs[alev],
offset);
1854template <
typename MF>
1858 auto const&
offset = linop.getSolvabilityOffset(amrlev, mglev, mf);
1860 for (
int c = 0; c < ncomp; ++c) {
1862 <<
" from mf component c = " << c
1863 <<
" on level (" << amrlev <<
", " << mglev <<
")\n";
1866 linop.fixSolvabilityByOffset(amrlev, mglev, mf,
offset);
1869#if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
1870template <
typename MF>
1871template <
class TMF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,
int>>
1875 const int amrlev = 0;
1876 const int mglev = linop.NMGLevels(amrlev) - 1;
1880 if (linop.isCellCentered())
1882 if (hypre_solver ==
nullptr)
1884 hypre_solver = linop.makeHypre(hypre_interface);
1886 hypre_solver->setVerbose(bottom_verbose);
1888 hypre_solver->setHypreOptionsNamespace(hypre_options_namespace);
1890 hypre_solver->setHypreOldDefault(hypre_old_default);
1891 hypre_solver->setHypreRelaxType(hypre_relax_type);
1892 hypre_solver->setHypreRelaxOrder(hypre_relax_order);
1893 hypre_solver->setHypreNumSweeps(hypre_num_sweeps);
1894 hypre_solver->setHypreStrongThreshold(hypre_strong_threshold);
1897 const BoxArray& ba = linop.m_grids[amrlev].back();
1898 const DistributionMapping& dm = linop.m_dmap[amrlev].back();
1899 const Geometry& geom = linop.m_geom[amrlev].back();
1901 hypre_bndry = std::make_unique<MLMGBndryT<MF>>(ba, dm, ncomp, geom);
1902 hypre_bndry->setHomogValues();
1903 const Real* dx = linop.m_geom[0][0].CellSize();
1904 IntVect crse_ratio = linop.m_coarse_data_crse_ratio.
allGT(0) ? linop.m_coarse_data_crse_ratio :
IntVect(1);
1905 RealVect bclocation(
AMREX_D_DECL(0.5*dx[0]*crse_ratio[0],
1906 0.5*dx[1]*crse_ratio[1],
1907 0.5*dx[2]*crse_ratio[2]));
1908 hypre_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc,
IntVect(-1), bclocation,
1909 linop.m_coarse_fine_bc_type);
1913 amrex::Real hypre_abstol =
1915 ? bottom_abstol : Real(-1.0);
1916 hypre_solver->solve(
1917 x,
b, bottom_reltol, hypre_abstol, bottom_maxiter, *hypre_bndry,
1918 linop.getMaxOrder());
1922 if (hypre_node_solver ==
nullptr)
1925 linop.makeHypreNodeLap(bottom_verbose, hypre_options_namespace);
1927 hypre_node_solver->solve(
x,
b, bottom_reltol, bottom_abstol, bottom_maxiter);
1932 if (linop.isSingular(amrlev) && linop.getEnforceSingularSolvable())
1934 makeSolvable(amrlev, mglev,
x);
1939#if defined(AMREX_USE_PETSC) && (AMREX_SPACEDIM > 1)
1940template <
typename MF>
1941template <
class TMF,std::enable_if_t<std::is_same_v<TMF,MultiFab>,
int>>
1943MLMGT<MF>::bottomSolveWithPETSc (MF&
x,
const MF&
b)
1947 if(petsc_solver ==
nullptr)
1949 petsc_solver = linop.makePETSc();
1950 petsc_solver->setVerbose(bottom_verbose);
1952 const BoxArray& ba = linop.m_grids[0].back();
1953 const DistributionMapping& dm = linop.m_dmap[0].back();
1954 const Geometry& geom = linop.m_geom[0].back();
1956 petsc_bndry = std::make_unique<MLMGBndryT<MF>>(ba, dm, ncomp, geom);
1957 petsc_bndry->setHomogValues();
1958 const Real* dx = linop.m_geom[0][0].CellSize();
1959 auto crse_ratio = linop.m_coarse_data_crse_ratio.allGT(0) ? linop.m_coarse_data_crse_ratio :
IntVect(1);
1960 RealVect bclocation(
AMREX_D_DECL(0.5*dx[0]*crse_ratio[0],
1961 0.5*dx[1]*crse_ratio[1],
1962 0.5*dx[2]*crse_ratio[2]));
1963 petsc_bndry->setLOBndryConds(linop.m_lobc, linop.m_hibc,
IntVect(-1), bclocation,
1964 linop.m_coarse_fine_bc_type);
1966 petsc_solver->solve(
x,
b, bottom_reltol, Real(-1.), bottom_maxiter, *petsc_bndry,
1967 linop.getMaxOrder());
1971template <
typename MF>
1975 RT a_tol_rel,
RT a_tol_abs,
const char* a_file_name)
const
1977 std::string file_name(a_file_name);
1982 std::string HeaderFileName(std::string(a_file_name)+
"/Header");
1983 std::ofstream HeaderFile;
1984 HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out |
1985 std::ofstream::trunc |
1986 std::ofstream::binary);
1987 if( ! HeaderFile.good()) {
1991 HeaderFile.precision(17);
1993 HeaderFile << linop.name() <<
"\n"
1994 <<
"a_tol_rel = " << a_tol_rel <<
"\n"
1995 <<
"a_tol_abs = " << a_tol_abs <<
"\n"
1996 <<
"verbose = " <<
verbose <<
"\n"
1997 <<
"max_iters = " << max_iters <<
"\n"
1998 <<
"nu1 = " << nu1 <<
"\n"
1999 <<
"nu2 = " << nu2 <<
"\n"
2000 <<
"nuf = " << nuf <<
"\n"
2001 <<
"nub = " << nub <<
"\n"
2002 <<
"max_fmg_iters = " << max_fmg_iters <<
"\n"
2003 <<
"bottom_solver = " <<
static_cast<int>(bottom_solver) <<
"\n"
2004 <<
"bottom_verbose = " << bottom_verbose <<
"\n"
2005 <<
"bottom_maxiter = " << bottom_maxiter <<
"\n"
2006 <<
"bottom_reltol = " << bottom_reltol <<
"\n"
2007 <<
"always_use_bnorm = " << always_use_bnorm <<
"\n"
2008 <<
"namrlevs = " << namrlevs <<
"\n"
2009 <<
"finest_amr_lev = " << finest_amr_lev <<
"\n"
2010 <<
"linop_prepared = " << linop_prepared <<
"\n"
2011 <<
"solve_called = " << solve_called <<
"\n";
2013 for (
int ilev = 0; ilev <= finest_amr_lev; ++ilev) {
2020 for (
int ilev = 0; ilev <= finest_amr_lev; ++ilev) {
2021 VisMF::Write(*a_sol[ilev], file_name+
"/Level_"+std::to_string(ilev)+
"/sol");
2022 VisMF::Write(*a_rhs[ilev], file_name+
"/Level_"+std::to_string(ilev)+
"/rhs");
2025 linop.checkPoint(file_name+
"/linop");
2028template <
typename MF>
2032 print_ident.resize(print_ident.size()+4,
' ');
2035template <
typename MF>
2039 if (print_ident.size() > 4) {
2040 print_ident.resize(print_ident.size()-4,
' ');
2042 print_ident.clear();
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define BL_PROFILE_VAR_STOP(vname)
Definition AMReX_BLProfiler.H:563
#define BL_PROFILE_VAR(fname, vname)
Definition AMReX_BLProfiler.H:560
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
Print on all processors of the default communicator.
Definition AMReX_Print.H:117
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:550
static bool SameRefs(const BoxArray &lhs, const BoxArray &rhs)
whether two BoxArrays share the same data
Definition AMReX_BoxArray.H:820
Calculates the distribution of FABs to MPI processes.
Definition AMReX_DistributionMapping.H:41
static bool SameRefs(const DistributionMapping &lhs, const DistributionMapping &rhs)
Definition AMReX_DistributionMapping.H:187
Definition AMReX_EBFabFactory.H:24
Solve using GMRES with multigrid as preconditioner.
Definition AMReX_GMRES_MLMG.H:21
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:73
Periodicity periodicity() const noexcept
Definition AMReX_Geometry.H:355
Interface
Definition AMReX_Hypre.H:21
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool allGT(const IntVectND< dim > &rhs) const noexcept
Returns true if this is greater than argument for all components. NOTE: This is NOT a strict weak ord...
Definition AMReX_IntVect.H:418
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr IntVectND< dim > TheDimensionVector(int d) noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:691
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE IntVectND< new_dim > resize(int fill_extra=0) const noexcept
Returns a new IntVectND of size new_dim by either shrinking or expanding this IntVectND.
Definition AMReX_IntVect.H:768
Definition AMReX_MLCGSolver.H:12
void setSolver(Type _typ) noexcept
Definition AMReX_MLCGSolver.H:28
void setVerbose(int _verbose)
Definition AMReX_MLCGSolver.H:39
int getNumIters() const noexcept
Definition AMReX_MLCGSolver.H:64
void setInitSolnZeroed(bool _sol_zeroed)
Definition AMReX_MLCGSolver.H:53
void setPrintIdentation(std::string s)
Definition AMReX_MLCGSolver.H:45
int solve(MF &solnL, const MF &rhsL, RT eps_rel, RT eps_abs)
Definition AMReX_MLCGSolver.H:89
Type
Definition AMReX_MLCGSolver.H:18
void setNGhost(int _nghost)
Definition AMReX_MLCGSolver.H:56
void setMaxIter(int _maxiter)
Definition AMReX_MLCGSolver.H:42
Definition AMReX_MLLinOp.H:98
typename FabDataType< MF >::fab_type FAB
Definition AMReX_MLLinOp.H:108
typename FabDataType< MF >::value_type RT
Definition AMReX_MLLinOp.H:109
Definition AMReX_MLMG.H:17
Definition AMReX_MLMG.H:12
void prepareForFluxes(Vector< MF const * > const &a_sol)
Definition AMReX_MLMG.H:579
void setBottomVerbose(int v) noexcept
Definition AMReX_MLMG.H:143
void setMaxFmgIter(int n) noexcept
Definition AMReX_MLMG.H:131
RT MLResNormInf(int alevmax, bool local=false)
Definition AMReX_MLMG.H:1811
RT MLRhsNormInf(bool local=false)
Definition AMReX_MLMG.H:1826
MLMGT(MLMGT< MF > &&)=delete
void actualBottomSolve()
Definition AMReX_MLMG.H:1495
Vector< Vector< MF > > rescor
Definition AMReX_MLMG.H:331
MF MFType
Definition AMReX_MLMG.H:25
BottomSolver getBottomSolver() const noexcept
Definition AMReX_MLMG.H:141
bool linop_prepared
Definition AMReX_MLMG.H:277
int bottom_verbose
Definition AMReX_MLMG.H:263
void checkPoint(const Vector< MultiFab * > &a_sol, const Vector< MultiFab const * > &a_rhs, RT a_tol_rel, RT a_tol_abs, const char *a_file_name) const
Definition AMReX_MLMG.H:1973
void setPreSmooth(int n) noexcept
Definition AMReX_MLMG.H:135
void setBottomToleranceAbs(RT t) noexcept
Definition AMReX_MLMG.H:146
Vector< Vector< MF > > cor
= rhs - L(sol)
Definition AMReX_MLMG.H:329
RT getFinalResidual() const noexcept
Definition AMReX_MLMG.H:236
int do_nsolve
N Solve.
Definition AMReX_MLMG.H:281
void interpCorrection(int alev)
Definition AMReX_MLMG.H:1713
Long solve_called
Definition AMReX_MLMG.H:278
void setBottomSmooth(int n) noexcept
Definition AMReX_MLMG.H:138
int final_fill_bc
Definition AMReX_MLMG.H:270
bool precond_mode
Definition AMReX_MLMG.H:246
void setNSolve(int flag) noexcept
Definition AMReX_MLMG.H:155
int getBottomVerbose() const
Definition AMReX_MLMG.H:123
std::unique_ptr< MF > ns_sol
Definition AMReX_MLMG.H:285
std::string print_ident
Definition AMReX_MLMG.H:288
Vector< Vector< MF > > cor_hold
Definition AMReX_MLMG.H:330
void computeResOfCorrection(int amrlev, int mglev)
Definition AMReX_MLMG.H:1790
void applyPrecond(const Vector< MF * > &out, const Vector< MF * > &in)
out = L(in) as a preconditioner
Definition AMReX_MLMG.H:970
void setCFStrategy(CFStrategy a_cf_strategy) noexcept
Definition AMReX_MLMG.H:142
void computeResWithCrseCorFineCor(int falev)
Definition AMReX_MLMG.H:1689
void NSolve(MLMGT< MF > &a_solver, MF &a_sol, MF &a_rhs)
Definition AMReX_MLMG.H:1471
typename MLLinOpT< MF >::Location Location
Definition AMReX_MLMG.H:30
RT m_final_resnorm0
Definition AMReX_MLMG.H:339
void apply(const Vector< MF * > &out, const Vector< MF * > &in)
out = L(in). Note that, if no actual solve is needed, one could turn off multigrid coarsening by cons...
Definition AMReX_MLMG.H:894
Vector< MF > rhs
Definition AMReX_MLMG.H:319
void setNSolveGridSize(int s) noexcept
Definition AMReX_MLMG.H:156
int nuf
when smoother is used as bottom solver
Definition AMReX_MLMG.H:256
void setVerbose(int v) noexcept
Definition AMReX_MLMG.H:129
int nu1
pre
Definition AMReX_MLMG.H:254
void computeMLResidual(int amrlevmax)
Definition AMReX_MLMG.H:1628
Vector< RT > m_iter_fine_resnorm0
Definition AMReX_MLMG.H:341
RT getInitResidual() const noexcept
Definition AMReX_MLMG.H:234
int getNumIters() const noexcept
Definition AMReX_MLMG.H:239
void setPostSmooth(int n) noexcept
Definition AMReX_MLMG.H:136
BottomSolver bottom_solver
Definition AMReX_MLMG.H:261
void mgVcycle(int amrlev, int mglev)
Definition AMReX_MLMG.H:1301
MLLinOpT< MF > & linop
Definition AMReX_MLMG.H:272
int ncomp
Definition AMReX_MLMG.H:273
void prepareForNSolve()
Definition AMReX_MLMG.H:1188
RT precond(Vector< MF * > const &a_sol, Vector< MF const * > const &a_rhs, RT a_tol_rel, RT a_tol_abs)
Definition AMReX_MLMG.H:561
int always_use_bnorm
Definition AMReX_MLMG.H:268
int max_fmg_iters
Definition AMReX_MLMG.H:259
void makeSolvable()
Definition AMReX_MLMG.H:1840
void setBottomSolver(BottomSolver s) noexcept
Definition AMReX_MLMG.H:140
int namrlevs
Definition AMReX_MLMG.H:274
void preparePrecond()
Definition AMReX_MLMG.H:1180
int nub
additional smoothing after bottom cg solver
Definition AMReX_MLMG.H:257
void incPrintIdentation()
Definition AMReX_MLMG.H:2030
typename MLLinOpT< MF >::RT RT
Definition AMReX_MLMG.H:27
int max_iters
Definition AMReX_MLMG.H:250
int do_fixed_number_of_iters
Definition AMReX_MLMG.H:251
void setThrowException(bool t) noexcept
Definition AMReX_MLMG.H:128
void decPrintIdentation()
Definition AMReX_MLMG.H:2037
void setFixedIter(int nit) noexcept
Definition AMReX_MLMG.H:132
int finest_amr_lev
Definition AMReX_MLMG.H:275
void getGradSolution(const Vector< Array< AMF *, AMREX_SPACEDIM > > &a_grad_sol, Location a_loc=Location::FaceCenter)
Definition AMReX_MLMG.H:590
Vector< RT > const & getResidualHistory() const noexcept
Definition AMReX_MLMG.H:238
void prepareLinOp()
Definition AMReX_MLMG.H:1168
RT m_rhsnorm0
Definition AMReX_MLMG.H:337
void setPrecondIter(int nit) noexcept
Definition AMReX_MLMG.H:133
CFStrategy
Definition AMReX_MLMG.H:33
void prepareForSolve(Vector< AMF * > const &a_sol, Vector< AMF const * > const &a_rhs)
Definition AMReX_MLMG.H:982
int bottomSolveWithCG(MF &x, const MF &b, typename MLCGSolverT< MF >::Type type)
Definition AMReX_MLMG.H:1607
std::unique_ptr< MLMGT< MF > > ns_mlmg
Definition AMReX_MLMG.H:284
void setAlwaysUseBNorm(int flag) noexcept
Definition AMReX_MLMG.H:149
void compResidual(const Vector< MF * > &a_res, const Vector< MF * > &a_sol, const Vector< MF const * > &a_rhs)
Definition AMReX_MLMG.H:826
void miniCycle(int amrlev)
Definition AMReX_MLMG.H:1290
RT solve(std::initializer_list< AMF * > a_sol, std::initializer_list< AMF const * > a_rhs, RT a_tol_rel, RT a_tol_abs, const char *checkpoint_file=nullptr)
int bottom_maxiter
Definition AMReX_MLMG.H:264
bool throw_exception
Definition AMReX_MLMG.H:247
RT bottom_reltol
Definition AMReX_MLMG.H:265
Vector< Vector< MF > > res
First Vector: Amr levels. 0 is the coarest level Second Vector: MG levels. 0 is the finest level.
Definition AMReX_MLMG.H:328
void setFinalFillBC(int flag) noexcept
Definition AMReX_MLMG.H:151
typename MLLinOpT< MF >::BCMode BCMode
Definition AMReX_MLMG.H:29
void computeResWithCrseSolFineCor(int calev, int falev)
Definition AMReX_MLMG.H:1656
MLMGT< MF > & operator=(MLMGT< MF > const &)=delete
Vector< MF > sol
Hypre.
Definition AMReX_MLMG.H:318
int nsolve_grid_size
Definition AMReX_MLMG.H:282
MLMGT(MLLinOpT< MF > &a_lp)
Definition AMReX_MLMG.H:350
int verbose
Definition AMReX_MLMG.H:248
void computeResidual(int alev)
Definition AMReX_MLMG.H:1646
MLLinOpT< MF > & getLinOp()
Definition AMReX_MLMG.H:242
std::unique_ptr< MF > ns_rhs
Definition AMReX_MLMG.H:286
timer_types
Definition AMReX_MLMG.H:334
@ ntimers
Definition AMReX_MLMG.H:334
@ bottom_time
Definition AMReX_MLMG.H:334
@ iter_time
Definition AMReX_MLMG.H:334
@ solve_time
Definition AMReX_MLMG.H:334
CFStrategy cf_strategy
Definition AMReX_MLMG.H:262
typename MLLinOpT< MF >::FAB FAB
Definition AMReX_MLMG.H:26
Vector< int > sol_is_alias
Definition AMReX_MLMG.H:322
RT getBottomToleranceAbs() const noexcept
Definition AMReX_MLMG.H:147
int numAMRLevels() const noexcept
Definition AMReX_MLMG.H:153
RT bottom_abstol
Definition AMReX_MLMG.H:266
Vector< int > m_niters_cg
Definition AMReX_MLMG.H:340
MLMGT(MLMGT< MF > const &)=delete
void mgFcycle()
Definition AMReX_MLMG.H:1415
Vector< double > timer
Definition AMReX_MLMG.H:335
RT getInitRHS() const noexcept
Definition AMReX_MLMG.H:232
RT ResNormInf(int alev, bool local=false)
Definition AMReX_MLMG.H:1802
Vector< int > const & getNumCGIters() const noexcept
Definition AMReX_MLMG.H:240
int nu2
post
Definition AMReX_MLMG.H:255
int max_precond_iters
Definition AMReX_MLMG.H:252
void bottomSolve()
Definition AMReX_MLMG.H:1457
void getFluxes(const Vector< Array< AMF *, AMREX_SPACEDIM > > &a_flux, Location a_loc=Location::FaceCenter)
For (alpha * a - beta * (del dot b grad)) phi = rhs, flux means -b grad phi
Definition AMReX_MLMG.H:621
void setBottomTolerance(RT t) noexcept
Definition AMReX_MLMG.H:145
void setFinalSmooth(int n) noexcept
Definition AMReX_MLMG.H:137
std::unique_ptr< MLLinOpT< MF > > ns_linop
Definition AMReX_MLMG.H:283
void addInterpCorrection(int alev, int mglev)
Definition AMReX_MLMG.H:1760
RT m_init_resnorm0
Definition AMReX_MLMG.H:338
int getVerbose() const
Definition AMReX_MLMG.H:122
RT solve(const Vector< AMF * > &a_sol, const Vector< AMF const * > &a_rhs, RT a_tol_rel, RT a_tol_abs, const char *checkpoint_file=nullptr)
void setBottomMaxIter(int n) noexcept
Definition AMReX_MLMG.H:144
void oneIter(int iter)
Definition AMReX_MLMG.H:1221
void setMaxIter(int n) noexcept
Definition AMReX_MLMG.H:130
This class provides the user with a few print options.
Definition AMReX_Print.H:35
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:27
Long size() const noexcept
Definition AMReX_Vector.H:50
static Long Write(const FabArray< FArrayBox > &mf, const std::string &name, VisMF::How how=NFiles, bool set_ghost=false)
Write a FabArray<FArrayBox> to disk in a "smart" way. Returns the total number of bytes written on th...
Definition AMReX_VisMF.cpp:975
void Max(KeyValuePair< K, V > &vi, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:126
void push(MPI_Comm c)
Definition AMReX_ParallelContext.H:102
void BarrierSub() noexcept
Definition AMReX_ParallelContext.H:88
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
int MyProcSub() noexcept
my sub-rank in current frame
Definition AMReX_ParallelContext.H:76
void pop()
Note that it's the user's responsibility to free the MPI_Comm.
Definition AMReX_ParallelContext.H:108
bool IOProcessorSub() noexcept
Am IO processor for current frame?
Definition AMReX_ParallelContext.H:80
Definition AMReX_Amr.cpp:49
MF::value_type norminf(MF const &mf, int scomp, int ncomp, IntVect const &nghost, bool local=false)
Definition AMReX_FabArrayUtility.H:1883
int nComp(FabArrayBase const &fa)
void FileOpenFailed(const std::string &file)
Output a message and abort when couldn't open the file.
Definition AMReX_Utility.cpp:131
DistributionMapping const & DistributionMap(FabArrayBase const &fa)
void average_down(const MultiFab &S_fine, MultiFab &S_crse, const Geometry &fgeom, const Geometry &cgeom, int scomp, int ncomp, int rr)
Definition AMReX_MultiFabUtil.cpp:314
IntVect nGrowVect(FabArrayBase const &fa)
void EB_average_face_to_cellcenter(MultiFab &ccmf, int dcomp, const Array< MultiFab const *, AMREX_SPACEDIM > &fmf)
Definition AMReX_EBMultiFabUtil.cpp:806
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > convert(const BoxND< dim > &b, const IntVectND< dim > &typ) noexcept
Returns a BoxND with different type.
Definition AMReX_Box.H:1435
void LocalCopy(DMF &dst, SMF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst = src
Definition AMReX_FabArrayUtility.H:1831
void EB_set_covered(MultiFab &mf, Real val)
Definition AMReX_EBMultiFabUtil.cpp:21
void LocalAdd(MF &dst, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &nghost)
dst += src
Definition AMReX_FabArrayUtility.H:1839
double second() noexcept
Definition AMReX_Utility.cpp:922
void UtilCreateCleanDirectory(const std::string &path, bool callbarrier=true)
Create a new directory, renaming the old one if it exists.
Definition AMReX_Utility.cpp:161
void EB_average_down(const MultiFab &S_fine, MultiFab &S_crse, const MultiFab &vol_fine, const MultiFab &vfrac_fine, int scomp, int ncomp, const IntVect &ratio)
Definition AMReX_EBMultiFabUtil.cpp:336
Vector< std::array< T *, AMREX_SPACEDIM > > GetVecOfArrOfPtrs(const Vector< std::array< std::unique_ptr< T >, AMREX_SPACEDIM > > &a)
Definition AMReX_Vector.H:138
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
BottomSolver
Definition AMReX_MLLinOp.H:30
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:111
void average_face_to_cellcenter(MultiFab &cc, int dcomp, const Vector< const MultiFab * > &fc, int ngrow)
Average face-based MultiFab onto cell-centered MultiFab.
Definition AMReX_MultiFabUtil.cpp:141
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T norm(const GpuComplex< T > &a_z) noexcept
Return the norm (magnitude squared) of a complex number.
Definition AMReX_GpuComplex.H:344
void setBndry(MF &dst, typename MF::value_type val, int scomp, int ncomp)
dst = val in ghost cells.
Definition AMReX_FabArrayUtility.H:1815
Vector< T * > GetVecOfPtrs(Vector< T > &a)
Definition AMReX_Vector.H:61
std::array< T *, AMREX_SPACEDIM > GetArrOfPtrs(std::array< T, AMREX_SPACEDIM > &a) noexcept
Definition AMReX_Array.H:852
void Scale(MF &dst, typename MF::value_type val, int scomp, int ncomp, int nghost)
dst *= val
Definition AMReX_FabArrayUtility.H:1822
std::array< T const *, AMREX_SPACEDIM > GetArrOfConstPtrs(const std::array< T, AMREX_SPACEDIM > &a) noexcept
Definition AMReX_Array.H:864
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:225
void ParallelCopy(MF &dst, MF const &src, int scomp, int dcomp, int ncomp, IntVect const &ng_src=IntVect(0), IntVect const &ng_dst=IntVect(0), Periodicity const &period=Periodicity::NonPeriodic())
dst = src w/ MPI communication
Definition AMReX_FabArrayUtility.H:1873
int verbose
Definition AMReX_DistributionMapping.cpp:36
void setVal(MF &dst, typename MF::value_type val)
dst = val
Definition AMReX_FabArrayUtility.H:1808
BoxArray const & boxArray(FabArrayBase const &fa)
std::array< T, N > Array
Definition AMReX_Array.H:24
BCMode
Definition AMReX_MLLinOp.H:85
Location
Definition AMReX_MLLinOp.H:87
FabArray memory allocation information.
Definition AMReX_FabArray.H:66