From 75571e2dcbf2417529c5ed8e24113580e8e1f3f1 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Wed, 20 Dec 2023 19:04:57 -0800 Subject: [PATCH] MLMG: Use free functions instead of MF member functions (#3681) Note that the use of unqualified functions (e.g., setVal instead of amrex::setVal) is intentional. With ADL, these calls in MLMG could work with user defined data. --- Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H | 72 ++++---- Src/LinearSolvers/MLMG/AMReX_MLLinOp.H | 83 ++++++--- Src/LinearSolvers/MLMG/AMReX_MLMG.H | 204 +++++++++++----------- 3 files changed, 199 insertions(+), 160 deletions(-) diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H index ff9198215f..3bfab3c9f6 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H @@ -12,8 +12,8 @@ class MLCGSolverT { public: - using FAB = typename MF::fab_type; - using RT = typename MF::value_type; + using FAB = typename MLLinOpT::FAB; + using RT = typename MLLinOpT::RT; enum struct Type { BiCGStab, CG }; @@ -99,12 +99,12 @@ MLCGSolverT::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) { BL_PROFILE("MLCGSolver::bicgstab"); - const int ncomp = sol.nComp(); + const int ncomp = nComp(sol); - MF p = Lp.make(amrlev, mglev, sol.nGrowVect()); - MF r = Lp.make(amrlev, mglev, sol.nGrowVect()); - p.setVal(RT(0.0)); // Make sure all entries are initialized to avoid errors - r.setVal(RT(0.0)); + MF p = Lp.make(amrlev, mglev, nGrowVect(sol)); + MF r = Lp.make(amrlev, mglev, nGrowVect(sol)); + setVal(p, RT(0.0)); // Make sure all entries are initialized to avoid errors + setVal(r, RT(0.0)); MF rh = Lp.make(amrlev, mglev, nghost); MF v = Lp.make(amrlev, mglev, nghost); @@ -114,19 +114,19 @@ MLCGSolverT::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) MF sorig; if ( initial_vec_zeroed ) { - r.LocalCopy(rhs,0,0,ncomp,nghost); + LocalCopy(r,rhs,0,0,ncomp,nghost); } else { sorig = Lp.make(amrlev, mglev, nghost); Lp.correctionResidual(amrlev, mglev, r, sol, rhs, MLLinOpT::BCMode::Homogeneous); - sorig.LocalCopy(sol,0,0,ncomp,nghost); - sol.setVal(RT(0.0)); + LocalCopy(sorig,sol,0,0,ncomp,nghost); + setVal(sol, RT(0.0)); } // Then normalize Lp.normalize(amrlev, mglev, r); - rh.LocalCopy (r ,0,0,ncomp,nghost); + LocalCopy(rh, r, 0,0,ncomp,nghost); RT rnorm = norm_inf(r); const RT rnorm0 = rnorm; @@ -159,13 +159,13 @@ MLCGSolverT::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) } if ( iter == 1 ) { - p.LocalCopy(r,0,0,ncomp,nghost); + LocalCopy(p,r,0,0,ncomp,nghost); } else { const RT beta = (rho/rho_1)*(alpha/omega); - MF::Saxpy(p, -omega, v, 0, 0, ncomp, nghost); // p += -omega*v - MF::Xpay(p, beta, r, 0, 0, ncomp, nghost); // p = r + beta*p + Saxpy(p, -omega, v, 0, 0, ncomp, nghost); // p += -omega*v + Xpay(p, beta, r, 0, 0, ncomp, nghost); // p = r + beta*p } Lp.apply(amrlev, mglev, v, p, MLLinOpT::BCMode::Homogeneous, MLLinOpT::StateMode::Correction); Lp.normalize(amrlev, mglev, v); @@ -179,8 +179,8 @@ MLCGSolverT::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) { ret = 2; break; } - MF::Saxpy(sol, alpha, p, 0, 0, ncomp, nghost); // sol += alpha * p - MF::Saxpy(r, -alpha, v, 0, 0, ncomp, nghost); // r += -alpha * v + Saxpy(sol, alpha, p, 0, 0, ncomp, nghost); // sol += alpha * p + Saxpy(r, -alpha, v, 0, 0, ncomp, nghost); // r += -alpha * v rnorm = norm_inf(r); rnorm = norm_inf(r); @@ -216,8 +216,8 @@ MLCGSolverT::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) { ret = 3; break; } - MF::Saxpy(sol, omega, r, 0, 0, ncomp, nghost); // sol += omega * r - MF::Saxpy(r, -omega, t, 0, 0, ncomp, nghost); // r += -omega * t + Saxpy(sol, omega, r, 0, 0, ncomp, nghost); // sol += omega * r + Saxpy(r, -omega, t, 0, 0, ncomp, nghost); // r += -omega * t rnorm = norm_inf(r); @@ -257,14 +257,14 @@ MLCGSolverT::solve_bicgstab (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) ) { if ( !initial_vec_zeroed ) { - sol.LocalAdd(sorig, 0, 0, ncomp, nghost); + LocalAdd(sol, sorig, 0, 0, ncomp, nghost); } } else { - sol.setVal(RT(0.0)); + setVal(sol, RT(0.0)); if ( !initial_vec_zeroed ) { - sol.LocalAdd(sorig, 0, 0, ncomp, nghost); + LocalAdd(sol, sorig, 0, 0, ncomp, nghost); } } @@ -277,10 +277,10 @@ MLCGSolverT::solve_cg (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) { BL_PROFILE("MLCGSolver::cg"); - const int ncomp = sol.nComp(); + const int ncomp = nComp(sol); - MF p = Lp.make(amrlev, mglev, sol.nGrowVect()); - p.setVal(RT(0.0)); + MF p = Lp.make(amrlev, mglev, nGrowVect(sol)); + setVal(p, RT(0.0)); MF r = Lp.make(amrlev, mglev, nghost); MF q = Lp.make(amrlev, mglev, nghost); @@ -288,14 +288,14 @@ MLCGSolverT::solve_cg (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) MF sorig; if ( initial_vec_zeroed ) { - r.LocalCopy(rhs,0,0,ncomp,nghost); + LocalCopy(r,rhs,0,0,ncomp,nghost); } else { sorig = Lp.make(amrlev, mglev, nghost); Lp.correctionResidual(amrlev, mglev, r, sol, rhs, MLLinOpT::BCMode::Homogeneous); - sorig.LocalCopy(sol,0,0,ncomp,nghost); - sol.setVal(RT(0.0)); + LocalCopy(sorig,sol,0,0,ncomp,nghost); + setVal(sol, RT(0.0)); } RT rnorm = norm_inf(r); @@ -330,12 +330,12 @@ MLCGSolverT::solve_cg (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) } if (iter == 1) { - p.LocalCopy(r,0,0,ncomp,nghost); + LocalCopy(p,r,0,0,ncomp,nghost); } else { RT beta = rho/rho_1; - MF::Xpay(p, beta, r, 0, 0, ncomp, nghost); // p = r + beta * p + Xpay(p, beta, r, 0, 0, ncomp, nghost); // p = r + beta * p } Lp.apply(amrlev, mglev, q, p, MLLinOpT::BCMode::Homogeneous, MLLinOpT::StateMode::Correction); @@ -357,8 +357,8 @@ MLCGSolverT::solve_cg (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) << " rho " << rho << " alpha " << alpha << '\n'; } - MF::Saxpy(sol, alpha, p, 0, 0, ncomp, nghost); // sol += alpha * p - MF::Saxpy(r, -alpha, q, 0, 0, ncomp, nghost); // r += -alpha * q + Saxpy(sol, alpha, p, 0, 0, ncomp, nghost); // sol += alpha * p + Saxpy(r, -alpha, q, 0, 0, ncomp, nghost); // r += -alpha * q rnorm = norm_inf(r); if ( verbose > 2 ) @@ -393,14 +393,14 @@ MLCGSolverT::solve_cg (MF& sol, const MF& rhs, RT eps_rel, RT eps_abs) if ( ( ret == 0 || ret == 8 ) && (rnorm < rnorm0) ) { if ( !initial_vec_zeroed ) { - sol.LocalAdd(sorig, 0, 0, ncomp, nghost); + LocalAdd(sol, sorig, 0, 0, ncomp, nghost); } } else { - sol.setVal(RT(0.0)); + setVal(sol, RT(0.0)); if ( !initial_vec_zeroed ) { - sol.LocalAdd(sorig, 0, 0, ncomp, nghost); + LocalAdd(sol, sorig, 0, 0, ncomp, nghost); } } @@ -422,8 +422,8 @@ template auto MLCGSolverT::norm_inf (const MF& res, bool local) -> RT { - int ncomp = res.nComp(); - RT result = res.norminf(0,ncomp,IntVect(0),true); + int ncomp = nComp(res); + RT result = norminf(res,0,ncomp,IntVect(0),true); if (!local) { BL_PROFILE("MLCGSolver::ParallelAllReduce"); ParallelAllReduce::Max(result, Lp.BottomCommunicator()); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H index b8aa71eebd..f0dca07f3a 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H @@ -85,6 +85,15 @@ struct LinOpEnumType enum struct Location { FaceCenter, FaceCentroid, CellCenter, CellCentroid }; }; +template struct LinOpData {}; +// +template +struct LinOpData > > +{ + using fab_type = typename T::fab_type; + using value_type = typename T::value_type; +}; + template class MLMGT; template class MLCGSolverT; template class MLPoissonT; @@ -100,8 +109,8 @@ public: template friend class MLPoissonT; template friend class MLABecLaplacianT; - using FAB = typename MF::fab_type; - using RT = typename MF::value_type; + using FAB = typename LinOpData::fab_type; + using RT = typename LinOpData::value_type; using BCType = LinOpBCType; using BCMode = LinOpEnumType::BCMode; @@ -1375,13 +1384,18 @@ template void MLLinOpT::make (Vector >& mf, IntVect const& ng) const { - mf.clear(); - mf.resize(m_num_amr_levels); - for (int alev = 0; alev < m_num_amr_levels; ++alev) { - mf[alev].resize(m_num_mg_levels[alev]); - for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) { - mf[alev][mlev] = make(alev, mlev, ng); + if constexpr (IsMultiFabLike_v) { + mf.clear(); + mf.resize(m_num_amr_levels); + for (int alev = 0; alev < m_num_amr_levels; ++alev) { + mf[alev].resize(m_num_mg_levels[alev]); + for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) { + mf[alev][mlev] = make(alev, mlev, ng); + } } + } else { + amrex::ignore_unused(mf, ng); + amrex::Abort("MLLinOpT::make: how did we get here?"); } } @@ -1389,39 +1403,62 @@ template MF MLLinOpT::make (int amrlev, int mglev, IntVect const& ng) const { - return MF(amrex::convert(m_grids[amrlev][mglev], m_ixtype), - m_dmap[amrlev][mglev], getNComp(), ng, MFInfo(), - *m_factory[amrlev][mglev]); + if constexpr (IsMultiFabLike_v) { + return MF(amrex::convert(m_grids[amrlev][mglev], m_ixtype), + m_dmap[amrlev][mglev], getNComp(), ng, MFInfo(), + *m_factory[amrlev][mglev]); + } else { + amrex::ignore_unused(amrlev, mglev, ng); + amrex::Abort("MLLinOpT::make: how did we get here?"); + return {}; + } } template MF MLLinOpT::makeAlias (MF const& mf) const { - return MF(mf, amrex::make_alias, 0, mf.nComp()); + if constexpr (IsMultiFabLike_v) { + return MF(mf, amrex::make_alias, 0, mf.nComp()); + } else { + amrex::ignore_unused(mf); + amrex::Abort("MLLinOpT::makeAlias: how did we get here?"); + return {}; + } } template MF MLLinOpT::makeCoarseMG (int amrlev, int mglev, IntVect const& ng) const { - BoxArray cba = m_grids[amrlev][mglev]; - IntVect ratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[mglev]; - cba.coarsen(ratio); - cba.convert(m_ixtype); - return MF(cba, m_dmap[amrlev][mglev], getNComp(), ng); - + if constexpr (IsMultiFabLike_v) { + BoxArray cba = m_grids[amrlev][mglev]; + IntVect ratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[mglev]; + cba.coarsen(ratio); + cba.convert(m_ixtype); + return MF(cba, m_dmap[amrlev][mglev], getNComp(), ng); + } else { + amrex::ignore_unused(amrlev, mglev, ng); + amrex::Abort("MLLinOpT::makeCoarseMG: how did we get here?"); + return {}; + } } template MF MLLinOpT::makeCoarseAmr (int famrlev, IntVect const& ng) const { - BoxArray cba = m_grids[famrlev][0]; - IntVect ratio(AMRRefRatio(famrlev-1)); - cba.coarsen(ratio); - cba.convert(m_ixtype); - return MF(cba, m_dmap[famrlev][0], getNComp(), ng); + if constexpr (IsMultiFabLike_v) { + BoxArray cba = m_grids[famrlev][0]; + IntVect ratio(AMRRefRatio(famrlev-1)); + cba.coarsen(ratio); + cba.convert(m_ixtype); + return MF(cba, m_dmap[famrlev][0], getNComp(), ng); + } else { + amrex::ignore_unused(famrlev, ng); + amrex::Abort("MLLinOpT::makeCoarseAmr: how did we get here?"); + return {}; + } } template diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H index 84adba7dfd..9bfc2f0007 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLMG.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.H @@ -21,8 +21,8 @@ public: template friend class MLCGSolverT; - using FAB = typename MF::fab_type; - using RT = typename MF::value_type; + using FAB = typename MLLinOpT::FAB; + using RT = typename MLLinOpT::RT; using BCMode = typename MLLinOpT::BCMode; using Location = typename MLLinOpT::Location; @@ -507,7 +507,7 @@ MLMGT::solve (const Vector& a_sol, const Vector& a_rhs, for (int alev = 0; alev < namrlevs; ++alev) { if (!sol_is_alias[alev]) { - a_sol[alev]->LocalCopy(sol[alev], 0, 0, ncomp, ng_back); + LocalCopy(*a_sol[alev], sol[alev], 0, 0, ncomp, ng_back); } } @@ -541,11 +541,11 @@ MLMGT::getGradSolution (const Vector >& a_grad_so Array grad_sol; for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { auto const& amf = *(a_grad_sol[alev][idim]); - grad_sol[idim].define(amf.boxArray(), amf.DistributionMap(), ncomp, 0); + grad_sol[idim].define(boxArray(amf), DistributionMap(amf), ncomp, 0); } linop.compGrad(alev, GetArrOfPtrs(grad_sol), sol[alev], a_loc); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - a_grad_sol[alev][idim]->LocalCopy(grad_sol[idim], 0, 0, ncomp, IntVect(0)); + LocalCopy(*a_grad_sol[alev][idim], grad_sol[idim], 0, 0, ncomp, IntVect(0)); } } } @@ -578,13 +578,13 @@ MLMGT::getFluxes (const Vector >& a_flux, for (int ilev = 0; ilev < namrlevs; ++ilev) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { auto const& amf = *(a_flux[ilev][idim]); - fluxes[ilev][idim].define(amf.boxArray(), amf.DistributionMap(), ncomp, 0); + fluxes[ilev][idim].define(boxArray(amf), DistributionMap(amf), ncomp, 0); } } getFluxes(GetVecOfArrOfPtrs(fluxes), GetVecOfPtrs(sol), a_loc); for (int ilev = 0; ilev < namrlevs; ++ilev) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - a_flux[ilev][idim]->LocalCopy(fluxes[ilev][idim], 0, 0, ncomp, IntVect(0)); + LocalCopy(*a_flux[ilev][idim], fluxes[ilev][idim], 0, 0, ncomp, IntVect(0)); } } } @@ -618,14 +618,14 @@ MLMGT::getFluxes (const Vector >& a_flux, for (int ilev = 0; ilev < namrlevs; ++ilev) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { auto const& amf = *(a_flux[ilev][idim]); - fluxes[ilev][idim].define(amf.boxArray(), amf.DistributionMap(), ncomp, 0); + fluxes[ilev][idim].define(boxArray(amf), DistributionMap(amf), ncomp, 0); } - sol[ilev].LocalCopy(*a_sol[ilev], 0, 0, ncomp, sol[ilev].nGrowVect()); + LocalCopy(sol[ilev], *a_sol[ilev], 0, 0, ncomp, nGrowVect(sol[ilev])); } linop.getFluxes(GetVecOfArrOfPtrs(fluxes), GetVecOfPtrs(sol), a_loc); for (int ilev = 0; ilev < namrlevs; ++ilev) { for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - a_flux[ilev][idim]->LocalCopy(fluxes[ilev][idim], 0, 0, ncomp, IntVect(0)); + LocalCopy(*a_flux[ilev][idim], fluxes[ilev][idim], 0, 0, ncomp, IntVect(0)); } } } @@ -653,11 +653,11 @@ MLMGT::getFluxes (const Vector & a_flux, Location a_loc) Vector fluxes(namrlevs); for (int ilev = 0; ilev < namrlevs; ++ilev) { auto const& amf = *a_flux[ilev]; - fluxes[ilev].define(amf.boxArray(), amf.DistributionMap(), ncomp, 0); + fluxes[ilev].define(boxArray(amf), DistributionMap(amf), ncomp, 0); } getFluxes(GetVecOfPtrs(fluxes), GetVecOfPtrs(sol), a_loc); for (int ilev = 0; ilev < namrlevs; ++ilev) { - a_flux[ilev]->LocalCopy(fluxes[ilev], 0, 0, ncomp, IntVect(0)); + LocalCopy(*a_flux[ilev], fluxes[ilev], 0, 0, ncomp, IntVect(0)); } } } @@ -676,11 +676,11 @@ void MLMGT::getFluxes (const Vector & a_flux, const Vector& a_sol, Location /*a_loc*/) { - AMREX_ASSERT(a_flux[0]->nComp() >= AMREX_SPACEDIM); + AMREX_ASSERT(nComp(*a_flux[0]) >= AMREX_SPACEDIM); if constexpr (! std::is_same()) { for (int alev = 0; alev < namrlevs; ++alev) { - sol[alev].LocalCopy(*a_sol[alev], 0, 0, ncomp, sol[alev].nGrowVect()); + LocalCopy(sol[alev], *a_sol[alev], 0, 0, ncomp, nGrowVect(sol[alev])); } } @@ -718,11 +718,11 @@ MLMGT::getFluxes (const Vector & a_flux, Vector fluxes(namrlevs); for (int ilev = 0; ilev < namrlevs; ++ilev) { auto const& amf = *a_flux[ilev]; - fluxes[ilev].define(amf.boxArray(), amf.DistributionMap(), ncomp, 0); + fluxes[ilev].define(boxArray(amf), DistributionMap(amf), ncomp, 0); } linop.getFluxes(GetVecOfPtrs(fluxes), GetVecOfPtrs(sol)); for (int ilev = 0; ilev < namrlevs; ++ilev) { - a_flux[ilev]->LocalCopy(fluxes[ilev], 0, 0, ncomp, IntVect(0)); + LocalCopy(*a_flux[ilev], fluxes[ilev], 0, 0, ncomp, IntVect(0)); } } } @@ -779,7 +779,7 @@ MLMGT::compResidual (const Vector& a_res, const Vector& a_sol, sol_is_alias.resize(namrlevs,true); for (int alev = 0; alev < namrlevs; ++alev) { - if (cf_strategy == CFStrategy::ghostnodes || a_sol[alev]->nGrowVect() == ng_sol) + if (cf_strategy == CFStrategy::ghostnodes || nGrowVect(*a_sol[alev]) == ng_sol) { sol[alev] = linop.makeAlias(*a_sol[alev]); sol_is_alias[alev] = true; @@ -790,7 +790,7 @@ MLMGT::compResidual (const Vector& a_res, const Vector& a_sol, { sol[alev] = linop.make(alev, 0, ng_sol); } - sol[alev].LocalCopy(*a_sol[alev], 0, 0, ncomp, IntVect(0)); + LocalCopy(sol[alev], *a_sol[alev], 0, 0, ncomp, IntVect(0)); } } @@ -808,9 +808,9 @@ MLMGT::compResidual (const Vector& a_res, const Vector& a_sol, const MF* prhs = a_rhs[alev]; #if (AMREX_SPACEDIM != 3) int nghost = (cf_strategy == CFStrategy::ghostnodes) ? linop.getNGrow(alev) : 0; - MF rhstmp(prhs->boxArray(), prhs->DistributionMap(), ncomp, nghost, + MF rhstmp(boxArray(*prhs), DistributionMap(*prhs), ncomp, nghost, MFInfo(), *linop.Factory(alev)); - rhstmp.LocalCopy(*prhs, 0, 0, ncomp, IntVect(nghost)); + LocalCopy(rhstmp, *prhs, 0, 0, ncomp, IntVect(nghost)); linop.applyMetricTerm(alev, 0, rhstmp); linop.unimposeNeumannBC(alev, rhstmp); linop.applyInhomogNeumannTerm(alev, rhstmp); @@ -822,9 +822,9 @@ MLMGT::compResidual (const Vector& a_res, const Vector& a_sol, *a_res[alev+1], sol[alev+1], *a_rhs[alev+1]); if (linop.isCellCentered()) { #ifdef AMREX_USE_EB - amrex::EB_average_down(*a_res[alev+1], *a_res[alev], 0, ncomp, amrrr[alev]); + EB_average_down(*a_res[alev+1], *a_res[alev], 0, ncomp, amrrr[alev]); #else - amrex::average_down(*a_res[alev+1], *a_res[alev], 0, ncomp, amrrr[alev]); + average_down(*a_res[alev+1], *a_res[alev], 0, ncomp, amrrr[alev]); #endif } } @@ -858,7 +858,7 @@ MLMGT::apply (const Vector& out, const Vector& a_in) nghost = linop.getNGrow(alev); in[alev] = a_in[alev]; } - else if (a_in[alev]->nGrowVect() == ng_sol) + else if (nGrowVect(*a_in[alev]) == ng_sol) { in[alev] = a_in[alev]; } @@ -866,18 +866,18 @@ MLMGT::apply (const Vector& out, const Vector& a_in) { IntVect ng = ng_sol; if (cf_strategy == CFStrategy::ghostnodes) { ng = IntVect(nghost); } - in_raii[alev].define(a_in[alev]->boxArray(), - a_in[alev]->DistributionMap(), - a_in[alev]->nComp(), ng, + in_raii[alev].define(boxArray (*a_in[alev]), + DistributionMap(*a_in[alev]), + nComp (*a_in[alev]), ng, MFInfo(), *linop.Factory(alev)); - in_raii[alev].LocalCopy(*a_in[alev], 0, 0, ncomp, IntVect(nghost)); + LocalCopy(in_raii[alev], *a_in[alev], 0, 0, ncomp, IntVect(nghost)); in[alev] = &(in_raii[alev]); } - rh[alev].define(a_in[alev]->boxArray(), - a_in[alev]->DistributionMap(), - a_in[alev]->nComp(), nghost, MFInfo(), + rh[alev].define(boxArray (*a_in[alev]), + DistributionMap(*a_in[alev]), + nComp (*a_in[alev]), nghost, MFInfo(), *linop.Factory(alev)); - rh[alev].setVal(RT(0.0)); + setVal(rh[alev], RT(0.0)); } if (!linop_prepared) { @@ -901,9 +901,9 @@ MLMGT::apply (const Vector& out, const Vector& a_in) *out[alev+1], *in[alev+1], rh[alev+1]); if (linop.isCellCentered()) { #ifdef AMREX_USE_EB - amrex::EB_average_down(*out[alev+1], *out[alev], 0, out[alev]->nComp(), amrrr[alev]); + EB_average_down(*out[alev+1], *out[alev], 0, nComp(*out[alev]), amrrr[alev]); #else - amrex::average_down(*out[alev+1], *out[alev], 0, out[alev]->nComp(), amrrr[alev]); + average_down(*out[alev+1], *out[alev], 0, nComp(*out[alev]), amrrr[alev]); #endif } } @@ -970,10 +970,10 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& } else { - if (a_sol[alev]->nGrowVect() == ng_sol) { + if (nGrowVect(*a_sol[alev]) == ng_sol) { if constexpr (std::is_same()) { sol[alev] = linop.makeAlias(*a_sol[alev]); - sol[alev].setBndry(RT(0.0), 0, ncomp); + setBndry(sol[alev], RT(0.0), 0, ncomp); sol_is_alias[alev] = true; } } @@ -981,8 +981,8 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& if (!solve_called) { sol[alev] = linop.make(alev, 0, ng_sol); } - sol[alev].LocalCopy(*a_sol[alev], 0, 0, ncomp, IntVect(0)); - sol[alev].setBndry(RT(0.0), 0, ncomp); + LocalCopy(sol[alev], *a_sol[alev], 0, 0, ncomp, IntVect(0)); + setBndry(sol[alev], RT(0.0), 0, ncomp); } } } @@ -994,7 +994,7 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& if (!solve_called) { rhs[alev] = linop.make(alev, 0, ng_rhs); } - rhs[alev].LocalCopy(*a_rhs[alev], 0, 0, ncomp, ng_rhs); + LocalCopy(rhs[alev], *a_rhs[alev], 0, 0, ncomp, ng_rhs); linop.applyMetricTerm(alev, 0, rhs[alev]); linop.unimposeNeumannBC(alev, rhs[alev]); linop.applyInhomogNeumannTerm(alev, rhs[alev]); @@ -1036,8 +1036,8 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& const int nmglevs = linop.NMGLevels(alev); for (int mglev = 0; mglev < nmglevs; ++mglev) { - res [alev][mglev].setVal(RT(0.0)); - rescor[alev][mglev].setVal(RT(0.0)); + setVal(res [alev][mglev], RT(0.0)); + setVal(rescor[alev][mglev], RT(0.0)); } } @@ -1054,7 +1054,7 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& if (cf_strategy == CFStrategy::ghostnodes) { _ng=IntVect(linop.getNGrow(alev,mglev)); } cor[alev][mglev] = linop.make(alev, mglev, _ng); } - cor[alev][mglev].setVal(RT(0.0)); + setVal(cor[alev][mglev], RT(0.0)); } } @@ -1070,7 +1070,7 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& if (cf_strategy == CFStrategy::ghostnodes) { _ng=IntVect(linop.getNGrow(alev,mglev)); } cor_hold[alev][mglev] = linop.make(alev, mglev, _ng); } - cor_hold[alev][mglev].setVal(RT(0.0)); + setVal(cor_hold[alev][mglev], RT(0.0)); } } for (int alev = 1; alev < finest_amr_lev; ++alev) @@ -1081,7 +1081,7 @@ MLMGT::prepareForSolve (Vector const& a_sol, Vector const& if (cf_strategy == CFStrategy::ghostnodes) { _ng=IntVect(linop.getNGrow(alev)); } cor_hold[alev][0] = linop.make(alev, 0, _ng); } - cor_hold[alev][0].setVal(RT(0.0)); + setVal(cor_hold[alev][0], RT(0.0)); } if (linop.m_parent // no embedded N-Solve @@ -1110,30 +1110,32 @@ template void MLMGT::prepareForNSolve () { - ns_linop = linop.makeNLinOp(nsolve_grid_size); + if constexpr (IsMultiFabLike_v) { + ns_linop = linop.makeNLinOp(nsolve_grid_size); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(); } - - const BoxArray& ba = (*ns_linop).m_grids[0][0]; - const DistributionMapping& dm =(*ns_linop).m_dmap[0][0]; - - int ng = 1; - if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; } - ns_sol = std::make_unique(ba, dm, ncomp, ng, MFInfo(), *(ns_linop->Factory(0,0))); - ng = 0; - if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; } - ns_rhs = std::make_unique(ba, dm, ncomp, ng, MFInfo(), *(ns_linop->Factory(0,0))); - ns_sol->setVal(RT(0.0)); - ns_rhs->setVal(RT(0.0)); - - ns_linop->setLevelBC(0, ns_sol.get()); - - ns_mlmg = std::make_unique>(*ns_linop); - ns_mlmg->setVerbose(0); - ns_mlmg->setFixedIter(1); - ns_mlmg->setMaxFmgIter(20); - ns_mlmg->setBottomSolver(BottomSolver::smoother); + int nghost = 0; + if (cf_strategy == CFStrategy::ghostnodes) { nghost = linop.getNGrow(); } + + const BoxArray& ba = (*ns_linop).m_grids[0][0]; + const DistributionMapping& dm =(*ns_linop).m_dmap[0][0]; + + int ng = 1; + if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; } + ns_sol = std::make_unique(ba, dm, ncomp, ng, MFInfo(), *(ns_linop->Factory(0,0))); + ng = 0; + if (cf_strategy == CFStrategy::ghostnodes) { ng = nghost; } + ns_rhs = std::make_unique(ba, dm, ncomp, ng, MFInfo(), *(ns_linop->Factory(0,0))); + setVal(*ns_sol, RT(0.0)); + setVal(*ns_rhs, RT(0.0)); + + ns_linop->setLevelBC(0, ns_sol.get()); + + ns_mlmg = std::make_unique>(*ns_linop); + ns_mlmg->setVerbose(0); + ns_mlmg->setFixedIter(1); + ns_mlmg->setMaxFmgIter(20); + ns_mlmg->setBottomSolver(BottomSolver::smoother); + } } // in : Residual (res) on the finest AMR level @@ -1149,7 +1151,7 @@ void MLMGT::oneIter (int iter) IntVect nghost(0); if (cf_strategy == CFStrategy::ghostnodes) { nghost = IntVect(linop.getNGrow(alev)); } - sol[alev].LocalAdd(cor[alev][0], 0, 0, ncomp, nghost); + LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost); // compute residual for the coarse AMR level computeResWithCrseSolFineCor(alev-1,alev); @@ -1175,7 +1177,7 @@ void MLMGT::oneIter (int iter) IntVect nghost(0); if (cf_strategy == CFStrategy::ghostnodes) { nghost = IntVect(linop.getNGrow(0)); } - sol[0].LocalAdd(cor[0][0], 0, 0, ncomp, nghost); + LocalAdd(sol[0], cor[0][0], 0, 0, ncomp, nghost); } for (int alev = 1; alev <= finest_amr_lev; ++alev) @@ -1185,10 +1187,10 @@ void MLMGT::oneIter (int iter) IntVect nghost(0); if (cf_strategy == CFStrategy::ghostnodes) { nghost = IntVect(linop.getNGrow(alev)); } - sol[alev].LocalAdd(cor[alev][0], 0, 0, ncomp, nghost); + LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost); if (alev != finest_amr_lev) { - cor_hold[alev][0].LocalAdd(cor[alev][0], 0, 0, ncomp, nghost); + LocalAdd(cor_hold[alev][0], cor[alev][0], 0, 0, ncomp, nghost); } // Update fine AMR level correction @@ -1196,10 +1198,10 @@ void MLMGT::oneIter (int iter) miniCycle(alev); - sol[alev].LocalAdd(cor[alev][0], 0, 0, ncomp, nghost); + LocalAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost); if (alev != finest_amr_lev) { - cor[alev][0].LocalAdd(cor_hold[alev][0], 0, 0, ncomp, nghost); + LocalAdd(cor[alev][0], cor_hold[alev][0], 0, 0, ncomp, nghost); } } @@ -1231,12 +1233,12 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { - RT norm = res[amrlev][mglev].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(res[amrlev][mglev],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " DN: Norm before smooth " << norm << "\n"; } - cor[amrlev][mglev].setVal(RT(0.0)); + setVal(cor[amrlev][mglev], RT(0.0)); bool skip_fillboundary = true; for (int i = 0; i < nu1; ++i) { linop.smooth(amrlev, mglev, cor[amrlev][mglev], res[amrlev][mglev], skip_fillboundary); @@ -1248,7 +1250,7 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { - RT norm = rescor[amrlev][mglev].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(rescor[amrlev][mglev],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " DN: Norm after smooth " << norm << "\n"; } @@ -1262,7 +1264,7 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) { if (verbose >= 4) { - RT norm = res[amrlev][mglev_bottom].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(res[amrlev][mglev_bottom],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " DN: Norm before bottom " << norm << "\n"; } @@ -1270,7 +1272,7 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { computeResOfCorrection(amrlev, mglev_bottom); - RT norm = rescor[amrlev][mglev_bottom].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(rescor[amrlev][mglev_bottom],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " UP: Norm after bottom " << norm << "\n"; } @@ -1279,11 +1281,11 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) { if (verbose >= 4) { - RT norm = res[amrlev][mglev_bottom].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(res[amrlev][mglev_bottom],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " Norm before smooth " << norm << "\n"; } - cor[amrlev][mglev_bottom].setVal(RT(0.0)); + setVal(cor[amrlev][mglev_bottom], RT(0.0)); bool skip_fillboundary = true; for (int i = 0; i < nu1; ++i) { linop.smooth(amrlev, mglev_bottom, cor[amrlev][mglev_bottom], @@ -1293,7 +1295,7 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { computeResOfCorrection(amrlev, mglev_bottom); - RT norm = rescor[amrlev][mglev_bottom].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(rescor[amrlev][mglev_bottom],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " Norm after smooth " << norm << "\n"; } @@ -1308,7 +1310,7 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { computeResOfCorrection(amrlev, mglev); - RT norm = rescor[amrlev][mglev].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(rescor[amrlev][mglev],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " UP: Norm before smooth " << norm << "\n"; } @@ -1321,7 +1323,7 @@ MLMGT::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { computeResOfCorrection(amrlev, mglev); - RT norm = rescor[amrlev][mglev].norminf(0,ncomp,IntVect(0)); + RT norm = norminf(rescor[amrlev][mglev],0,ncomp,IntVect(0)); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " UP: Norm after smooth " << norm << "\n"; } @@ -1361,12 +1363,12 @@ MLMGT::mgFcycle () // rescor = res - L(cor) computeResOfCorrection(amrlev, mglev); // res = rescor; this provides b to the vcycle below - res[amrlev][mglev].LocalCopy(rescor[amrlev][mglev], 0, 0, ncomp, nghost); + LocalCopy(res[amrlev][mglev], rescor[amrlev][mglev], 0, 0, ncomp, nghost); // save cor; do v-cycle; add the saved to cor std::swap(cor[amrlev][mglev], cor_hold[amrlev][mglev]); mgVcycle(amrlev, mglev); - cor[amrlev][mglev].LocalAdd(cor_hold[amrlev][mglev], 0, 0, ncomp, nghost); + LocalAdd(cor[amrlev][mglev], cor_hold[amrlev][mglev], 0, 0, ncomp, nghost); } } @@ -1393,16 +1395,16 @@ MLMGT::NSolve (MLMGT& a_solver, MF& a_sol, MF& a_rhs) { BL_PROFILE("MLMG::NSolve()"); - a_sol.setVal(RT(0.0)); + setVal(a_sol, RT(0.0)); MF const& res_bottom = res[0].back(); - if (BoxArray::SameRefs(a_rhs.boxArray(),res_bottom.boxArray()) && - DistributionMapping::SameRefs(a_rhs.DistributionMap(),res_bottom.DistributionMap())) + if (BoxArray::SameRefs(boxArray(a_rhs),boxArray(res_bottom)) && + DistributionMapping::SameRefs(DistributionMap(a_rhs),DistributionMap(res_bottom))) { - a_rhs.LocalCopy(res_bottom, 0, 0, ncomp, IntVect(0)); + LocalCopy(a_rhs, res_bottom, 0, 0, ncomp, IntVect(0)); } else { - a_rhs.setVal(RT(0.0)); - a_rhs.ParallelCopy(res_bottom); + setVal(a_rhs, RT(0.0)); + ParallelCopy(a_rhs, res_bottom, 0, 0, ncomp); } a_solver.solve(Vector{&a_sol}, Vector{&a_rhs}, @@ -1428,7 +1430,7 @@ MLMGT::actualBottomSolve () auto& x = cor[amrlev][mglev]; auto& b = res[amrlev][mglev]; - x.setVal(RT(0.0)); + setVal(x, RT(0.0)); if (bottom_solver == BottomSolver::smoother) { @@ -1444,9 +1446,9 @@ MLMGT::actualBottomSolve () MF raii_b; if (linop.isBottomSingular() && linop.getEnforceSingularSolvable()) { - const IntVect ng = b.nGrowVect(); + const IntVect ng = nGrowVect(b); raii_b = linop.make(amrlev, mglev, ng); - raii_b.LocalCopy(b, 0, 0, ncomp, ng); + LocalCopy(raii_b, b, 0, 0, ncomp, ng); bottom_b = &raii_b; makeSolvable(amrlev,mglev,*bottom_b); @@ -1486,7 +1488,7 @@ MLMGT::actualBottomSolve () int ret = bottomSolveWithCG(x, *bottom_b, cg_type); // If the MLMG solve failed then set the correction to zero if (ret != 0) { - cor[amrlev][mglev].setVal(RT(0.0)); + setVal(cor[amrlev][mglev], RT(0.0)); if (bottom_solver == BottomSolver::cgbicg || bottom_solver == BottomSolver::bicgcg) { if (bottom_solver == BottomSolver::cgbicg) { @@ -1496,7 +1498,7 @@ MLMGT::actualBottomSolve () } ret = bottomSolveWithCG(x, *bottom_b, cg_type); if (ret != 0) { - cor[amrlev][mglev].setVal(RT(0.0)); + setVal(cor[amrlev][mglev], RT(0.0)); } else { // switch permanently if (cg_type == MLCGSolverT::Type::CG) { bottom_solver = BottomSolver::cg; @@ -1591,7 +1593,7 @@ MLMGT::computeResWithCrseSolFineCor (int calev, int falev) linop.solutionResidual(calev, crse_res, crse_sol, crse_rhs, crse_bcdata); linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res, BCMode::Homogeneous); - fine_res.LocalCopy(fine_rescor, 0, 0, ncomp, nghost); + LocalCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost); linop.reflux(calev, crse_res, crse_sol, crse_rhs, fine_res, fine_sol, fine_rhs); @@ -1619,7 +1621,7 @@ MLMGT::computeResWithCrseCorFineCor (int falev) // fine_rescor = fine_res - L(fine_cor) linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res, BCMode::Inhomogeneous, &crse_cor); - fine_res.LocalCopy(fine_rescor, 0, 0, ncomp, nghost); + LocalCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost); } // Interpolate correction from coarse to fine AMR level. @@ -1648,9 +1650,9 @@ MLMGT::interpCorrection (int alev) } MF cfine = linop.makeCoarseAmr(alev, IntVect(ng_dst)); - cfine.setVal(RT(0.0)); - cfine.ParallelCopy(crse_cor, 0, 0, ncomp, IntVect(ng_src), IntVect(ng_dst), - crse_geom.periodicity()); + setVal(cfine, RT(0.0)); + ParallelCopy(cfine, crse_cor, 0, 0, ncomp, IntVect(ng_src), IntVect(ng_dst), + crse_geom.periodicity()); linop.interpolationAmr(alev, fine_cor, cfine, nghost); // NOLINT(readability-suspicious-call-argument) } @@ -1689,7 +1691,7 @@ MLMGT::addInterpCorrection (int alev, int mglev) else { cfine = linop.makeCoarseMG(alev, mglev, IntVect(0)); - cfine.ParallelCopy(crse_cor,0,0,ncomp,IntVect(0),IntVect(0)); + ParallelCopy(cfine, crse_cor, 0, 0, ncomp); cmf = &cfine; }