From 3b50b5c0628571999d29720de15c1c2a20b97c54 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Thu, 2 Jun 2022 12:46:37 -0700 Subject: [PATCH] MLMG interface These changes are made to support a generic type (i.e., amrex::Any) in MLMG. This is still work in progress. But it should not break any existing codes. --- Src/Base/AMReX_Any.H | 5 +- .../MLMG/AMReX_MLABecLaplacian.cpp | 4 +- Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H | 5 + Src/LinearSolvers/MLMG/AMReX_MLCGSolver.cpp | 7 + Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H | 4 +- .../MLMG/AMReX_MLCellABecLap.cpp | 11 +- Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H | 27 +- Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp | 439 +++++++- .../MLMG/AMReX_MLEBNodeFDLaplacian.H | 2 +- .../MLMG/AMReX_MLEBNodeFDLaplacian.cpp | 21 +- Src/LinearSolvers/MLMG/AMReX_MLLinOp.H | 79 +- Src/LinearSolvers/MLMG/AMReX_MLLinOp.cpp | 290 +++++- Src/LinearSolvers/MLMG/AMReX_MLMG.H | 42 +- Src/LinearSolvers/MLMG/AMReX_MLMG.cpp | 962 +++++------------- .../MLMG/AMReX_MLNodeLaplacian.H | 8 +- .../MLMG/AMReX_MLNodeLaplacian.cpp | 29 +- Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H | 29 +- Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.cpp | 156 ++- 18 files changed, 1300 insertions(+), 820 deletions(-) diff --git a/Src/Base/AMReX_Any.H b/Src/Base/AMReX_Any.H index b57aa9a39ef..31c824825a4 100644 --- a/Src/Base/AMReX_Any.H +++ b/Src/Base/AMReX_Any.H @@ -60,15 +60,18 @@ public: private: struct innards_base { virtual const std::type_info& Type () const = 0; + virtual ~innards_base () = default; }; template struct innards : innards_base { - innards(MF && mf) + innards (MF && mf) : m_mf(std::forward(mf)) {} + virtual ~innards () = default; + virtual const std::type_info& Type () const override { return typeid(MF); } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.cpp b/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.cpp index 89dbb268e10..e5a9b0b31af 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.cpp @@ -323,10 +323,10 @@ MLABecLaplacian::applyMetricTermsCoeffs () for (int alev = 0; alev < m_num_amr_levels; ++alev) { const int mglev = 0; - applyMetricTerm(alev, mglev, m_a_coeffs[alev][mglev]); + applyMetricTermToMF(alev, mglev, m_a_coeffs[alev][mglev]); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - applyMetricTerm(alev, mglev, m_b_coeffs[alev][mglev][idim]); + applyMetricTermToMF(alev, mglev, m_b_coeffs[alev][mglev][idim]); } } #endif diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H index 45464bbeb9c..a33d70b4771 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.H @@ -40,6 +40,11 @@ public: Real eps_rel, Real eps_abs); + int solve (Any& solnL, + const Any& rhsL, + Real eps_rel, + Real eps_abs); + void setVerbose (int _verbose) { verbose = _verbose; } int getVerbose () const { return verbose; } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.cpp index c32b0d6199d..76144e6d42f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLCGSolver.cpp @@ -78,6 +78,13 @@ MLCGSolver::solve (MultiFab& sol, } } +int +MLCGSolver::solve (Any& sol, const Any& rhs, Real eps_rel, Real eps_abs) +{ + AMREX_ASSERT(sol.is()); // xxxxx TODO: MLCGSolver Any + return solve(sol.get(), rhs.get(), eps_rel, eps_abs); +} + int MLCGSolver::solve_bicgstab (MultiFab& sol, const MultiFab& rhs, diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H index 985bc9855b4..8849a2be292 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.H @@ -59,9 +59,9 @@ public: virtual MultiFab const* getACoeffs (int amrlev, int mglev) const = 0; virtual Array getBCoeffs (int amrlev, int mglev) const = 0; - virtual void applyInhomogNeumannTerm (int amrlev, MultiFab& rhs) const final override; + virtual void applyInhomogNeumannTerm (int amrlev, Any& rhs) const final override; - virtual void applyOverset (int amlev, MultiFab& rhs) const override; + virtual void applyOverset (int amlev, Any& rhs) const override; #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) virtual std::unique_ptr makeHypre (Hypre::Interface hypre_interface) const override; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp index b5580b3c15c..af094d89406 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellABecLap.cpp @@ -108,7 +108,7 @@ MLCellABecLap::define (const Vector& a_geom, amrlev = 0; for (int mglev = 1; mglev < m_num_mg_levels[amrlev]; ++mglev) { MultiFab foo(m_grids[amrlev][mglev], m_dmap[amrlev][mglev], 1, 0, MFInfo().SetAlloc(false)); - if (! isMFIterSafe(*m_overset_mask[amrlev][mglev], foo)) { + if (! amrex::isMFIterSafe(*m_overset_mask[amrlev][mglev], foo)) { auto osm = std::make_unique(m_grids[amrlev][mglev], m_dmap[amrlev][mglev], 1, 1); osm->ParallelCopy(*m_overset_mask[amrlev][mglev]); @@ -193,13 +193,16 @@ MLCellABecLap::getFluxes (const Vector >& a_flux } void -MLCellABecLap::applyInhomogNeumannTerm (int amrlev, MultiFab& rhs) const +MLCellABecLap::applyInhomogNeumannTerm (int amrlev, Any& a_rhs) const { bool has_inhomog_neumann = hasInhomogNeumannBC(); bool has_robin = hasRobinBC(); if (!has_inhomog_neumann && !has_robin) return; + AMREX_ASSERT(a_rhs.is()); + MultiFab& rhs = a_rhs.get(); + int ncomp = getNComp(); const int mglev = 0; @@ -414,9 +417,11 @@ MLCellABecLap::applyInhomogNeumannTerm (int amrlev, MultiFab& rhs) const } void -MLCellABecLap::applyOverset (int amrlev, MultiFab& rhs) const +MLCellABecLap::applyOverset (int amrlev, Any& a_rhs) const { if (m_overset_mask[amrlev][0]) { + AMREX_ASSERT(a_rhs.is()); + auto& rhs = a_rhs.get(); const int ncomp = getNComp(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H index f1168e5c41e..64d93851260 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H @@ -3,6 +3,7 @@ #include #include +#include namespace amrex { @@ -109,6 +110,8 @@ public: virtual void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const override; + virtual void interpAssign (int amrlev, int fmglev, MultiFab& fine, MultiFab& crse) const override; + virtual void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, const MultiFab& fine_sol, const MultiFab& fine_rhs) override; @@ -132,10 +135,15 @@ public: virtual void compGrad (int amrlev, const Array& grad, MultiFab& sol, Location loc) const override; - virtual void applyMetricTerm (int amrlev, int mglev, MultiFab& rhs) const final override; + virtual void applyMetricTerm (int amrlev, int mglev, Any& rhs) const final override; virtual void unapplyMetricTerm (int amrlev, int mglev, MultiFab& rhs) const final override; virtual void fillSolutionBC (int amrlev, MultiFab& sol, const MultiFab* crse_bcdata=nullptr) final override; + virtual Vector getSolvabilityOffset (int amrlev, int mglev, + Any const& rhs) const override; + virtual void fixSolvabilityByOffset (int amrlev, int mglev, Any& rhs, + Vector const& offset) const override; + virtual void prepareForSolve () override; virtual Real xdoty (int amrlev, int mglev, const MultiFab& x, const MultiFab& y, bool local) const final override; @@ -146,6 +154,18 @@ public: const Array& flux, const FArrayBox& sol, Location loc, const int face_only=0) const = 0; + // This could be turned into template if needed. + void applyMetricTermToMF (int amrlev, int mglev, MultiFab& rhs) const; + + virtual Real AnyNormInfMask (int amrlev, Any const& a, bool local) const override; + + virtual void AnyAvgDownResAmr (int clev, Any& cres, Any const& fres) const override; + + virtual void AnyInterpolationAmr (int famrlev, Any& fine, const Any& crse, + IntVect const& /*nghost*/) const override; + + virtual void AnyAverageDownAndSync (Vector& sol) const override; + struct BCTL { BoundCond type; Real location; @@ -210,12 +230,17 @@ protected: // boundary cell flags for covered, not_covered, outside_domain Vector > > m_maskvals; + Vector > m_norm_fine_mask; + mutable Vector m_fluxreg; private: void defineAuxData (); void defineBC (); + + void computeVolInv () const; + mutable Vector > m_volinv; // used by solvability fix }; } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp index 8f6921950e7..e7325972fd2 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.cpp @@ -1,6 +1,7 @@ #include #include +#include #include #ifndef BL_NO_FORT @@ -9,6 +10,11 @@ namespace amrex { +#ifdef AMREX_SOFT_PERF_COUNTERS +// perf_counters +MLCellLinOp::Counters MLCellLinOp::perf_counters; +#endif + namespace { // Have to put it here due to CUDA extended lambda limitation struct ABCTag { @@ -97,6 +103,7 @@ MLCellLinOp::defineAuxData () m_undrrelxr.resize(m_num_amr_levels); m_maskvals.resize(m_num_amr_levels); m_fluxreg.resize(m_num_amr_levels-1); + m_norm_fine_mask.resize(m_num_amr_levels-1); const int ncomp = getNComp(); @@ -136,6 +143,9 @@ MLCellLinOp::defineAuxData () m_dmap[amrlev+1][0], m_dmap[amrlev][0], m_geom[amrlev+1][0], m_geom[amrlev][0], ratio, amrlev+1, ncomp); + m_norm_fine_mask[amrlev] = std::make_unique + (makeFineMask(m_grids[amrlev][0], m_dmap[amrlev][0], m_grids[amrlev+1][0], + ratio, 1, 0)); } #if (AMREX_SPACEDIM != 3) @@ -1316,7 +1326,20 @@ MLCellLinOp::BndryCondLoc::setLOBndryConds (const Geometry& geom, const Real* dx } void -MLCellLinOp::applyMetricTerm (int amrlev, int mglev, MultiFab& rhs) const +MLCellLinOp::applyMetricTerm (int amrlev, int mglev, Any& rhs) const +{ + amrex::ignore_unused(amrlev,mglev,rhs); +#if (AMREX_SPACEDIM != 3) + + if (!m_has_metric_term) return; + + AMREX_ASSERT(rhs.is()); + applyMetricTermToMF(amrlev, mglev, rhs.get()); +#endif +} + +void +MLCellLinOp::applyMetricTermToMF (int amrlev, int mglev, MultiFab& rhs) const { amrex::ignore_unused(amrlev,mglev,rhs); #if (AMREX_SPACEDIM != 3) @@ -1435,9 +1458,417 @@ MLCellLinOp::update () if (MLLinOp::needsUpdate()) MLLinOp::update(); } -#ifdef AMREX_SOFT_PERF_COUNTERS -// perf_counters -MLCellLinOp::Counters MLCellLinOp::perf_counters; +void +MLCellLinOp::computeVolInv () const +{ + if (!m_volinv.empty()) return; + + m_volinv.resize(m_num_amr_levels); + for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev) { + m_volinv[amrlev].resize(NMGLevels(amrlev)); + } + + // We don't need to compute for every level + + auto f = [&] (int amrlev, int mglev) { +#ifdef AMREX_USE_EB + auto factory = dynamic_cast(Factory(amrlev,mglev)); + if (factory) + { + const MultiFab& vfrac = factory->getVolFrac(); + m_volinv[amrlev][mglev] = vfrac.sum(0,true); + } + else #endif + { + m_volinv[amrlev][mglev] + = Real(1.0 / compactify(Geom(amrlev,mglev).Domain()).d_numPts()); + } + }; + + // amrlev = 0, mglev = 0 + f(0,0); + + int mgbottom = NMGLevels(0)-1; + f(0,mgbottom); + +#ifdef AMREX_USE_EB + Real temp1, temp2; + auto factory = dynamic_cast(Factory(0,0)); + if (factory) + { + ParallelAllReduce::Sum({m_volinv[0][0], m_volinv[0][mgbottom]}, + ParallelContext::CommunicatorSub()); + temp1 = Real(1.0)/m_volinv[0][0]; + temp2 = Real(1.0)/m_volinv[0][mgbottom]; + } + else + { + temp1 = m_volinv[0][0]; + temp2 = m_volinv[0][mgbottom]; + } + m_volinv[0][0] = temp1; + m_volinv[0][mgbottom] = temp2; +#endif +} + +Vector +MLCellLinOp::getSolvabilityOffset (int amrlev, int mglev, Any const& a_rhs) const +{ + AMREX_ASSERT(a_rhs.is()); + auto const& rhs = a_rhs.get(); + + computeVolInv(); + + const int ncomp = getNComp(); + Vector offset(ncomp); + +#ifdef AMREX_USE_EB + auto factory = dynamic_cast(Factory(amrlev,mglev)); + if (factory) + { + const MultiFab& vfrac = factory->getVolFrac(); + for (int c = 0; c < ncomp; ++c) { + offset[c] = MultiFab::Dot(rhs, c, vfrac, 0, 1, 0, true) * m_volinv[amrlev][mglev]; + } + } + else +#endif + { + for (int c = 0; c < ncomp; ++c) { + offset[c] = rhs.sum(c,true) * m_volinv[amrlev][mglev]; + } + } + + ParallelAllReduce::Sum(offset.data(), ncomp, ParallelContext::CommunicatorSub()); + + return offset; +} + +Real +MLCellLinOp::AnyNormInfMask (int amrlev, Any const& a, bool local) const +{ + AMREX_ASSERT(a.is()); + auto& mf = a.get(); + + const int finest_level = NAMRLevels() - 1; + Real norm = 0._rt; +#ifdef AMREX_USE_EB + const int ncomp = getNComp(); + if (! mf.isAllRegular()) { + auto factory = dynamic_cast(Factory(amrlev)); + const MultiFab& vfrac = factory->getVolFrac(); + if (amrlev == finest_level) { +#ifdef AMREX_USE_GPU + if (Gpu::inLaunchRegion()) { + auto const& ma = mf.const_arrays(); + auto const& vfrac_ma = vfrac.const_arrays(); + norm = ParReduce(TypeList{}, TypeList{}, + mf, IntVect(0), ncomp, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) + -> GpuTuple + { + return amrex::Math::abs(ma[box_no](i,j,k,n) + *vfrac_ma[box_no](i,j,k)); + }); + } else +#endif + { +#ifdef AMREX_USE_OMP +#pragma omp parallel reduction(max:norm) +#endif + for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) { + Box const& bx = mfi.tilebox(); + auto const& fab = mf.const_array(mfi); + auto const& v = vfrac.const_array(mfi); + AMREX_LOOP_4D(bx, ncomp, i, j, k, n, + { + norm = std::max(norm, amrex::Math::abs(fab(i,j,k,n)*v(i,j,k))); + }); + } + } + } else { +#ifdef AMREX_USE_GPU + if (Gpu::inLaunchRegion()) { + auto const& ma = mf.const_arrays(); + auto const& mask_ma = m_norm_fine_mask[amrlev]->const_arrays(); + auto const& vfrac_ma = vfrac.const_arrays(); + norm = ParReduce(TypeList{}, TypeList{}, + mf, IntVect(0), ncomp, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) + -> GpuTuple + { + if (mask_ma[box_no](i,j,k)) { + return amrex::Math::abs(ma[box_no](i,j,k,n) + *vfrac_ma[box_no](i,j,k)); + } else { + return Real(0.0); + } + }); + } else +#endif + { +#ifdef AMREX_USE_OMP +#pragma omp parallel reduction(max:norm) +#endif + for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) { + Box const& bx = mfi.tilebox(); + auto const& fab = mf.const_array(mfi); + auto const& mask = m_norm_fine_mask[amrlev]->const_array(mfi); + auto const& v = vfrac.const_array(mfi); + AMREX_LOOP_4D(bx, ncomp, i, j, k, n, + { + if (mask(i,j,k)) { + norm = std::max(norm, amrex::Math::abs(fab(i,j,k,n)*v(i,j,k))); + } + }); + } + } + } + } else +#endif + { + iMultiFab const* fine_mask = (amrlev == finest_level) + ? nullptr : m_norm_fine_mask[amrlev].get(); + norm = MFNormInf(mf, fine_mask, true); + } + + if (!local) ParallelAllReduce::Max(norm, ParallelContext::CommunicatorSub()); + return norm; +} + +void +MLCellLinOp::AnyAvgDownResAmr (int clev, Any& cres, Any const& fres) const +{ + AMREX_ASSERT(cres.is() && fres.is()); +#ifdef AMREX_USE_EB + amrex::EB_average_down +#else + amrex::average_down +#endif + (fres.get(), cres.get(), 0, getNComp(), AMRRefRatio(clev)); +} + +void +MLCellLinOp::AnyInterpolationAmr (int famrlev, Any& a_fine, const Any& a_crse, + IntVect const& /*nghost*/) const +{ + AMREX_ASSERT(a_fine.is()); + MultiFab& fine = a_fine.get(); + MultiFab const& crse = a_crse.get(); + + const int ncomp = getNComp(); + const int refratio = AMRRefRatio(famrlev-1); + +#ifdef AMREX_USE_EB + auto factory = dynamic_cast(Factory(famrlev)); + const FabArray* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr; +#endif + + MFItInfo mfi_info; + if (Gpu::notInLaunchRegion()) mfi_info.EnableTiling().SetDynamic(true); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(fine, mfi_info); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + Array4 const& ff = fine.array(mfi); + Array4 const& cc = crse.const_array(mfi); +#ifdef AMREX_USE_EB + bool call_lincc; + if (factory) + { + const auto& flag = (*flags)[mfi]; + if (flag.getType(amrex::grow(bx,1)) == FabType::regular) { + call_lincc = true; + } else { + Array4 const& flg = flag.const_array(); + switch(refratio) { + case 2: + { + AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, + { + mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp); + }); + break; + } + case 4: + { + AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, + { + mlmg_eb_cc_interp_r<4>(tbx, ff, cc, flg, ncomp); + }); + break; + } + default: + amrex::Abort("mlmg_eb_cc_interp: only refratio 2 and 4 are supported"); + } + + call_lincc = false; + } + } + else + { + call_lincc = true; + } +#else + const bool call_lincc = true; +#endif + if (call_lincc) + { + switch(refratio) { + case 2: + { + AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, + { + mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp); + }); + break; + } + case 4: + { + AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, + { + mlmg_lin_cc_interp_r4(tbx, ff, cc, ncomp); + }); + break; + } + default: + amrex::Abort("mlmg_lin_cc_interp: only refratio 2 and 4 are supported"); + } + } + } +} + +void +MLCellLinOp::interpAssign (int amrlev, int fmglev, MultiFab& fine, MultiFab& crse) const +{ + const int ncomp = getNComp(); + + const Geometry& crse_geom = Geom(amrlev,fmglev+1); + const IntVect refratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[fmglev]; + const IntVect ng = crse.nGrowVect(); + + MultiFab cfine; + const MultiFab* cmf; + + if (amrex::isMFIterSafe(crse, fine)) + { + crse.FillBoundary(crse_geom.periodicity()); + cmf = &crse; + } + else + { + BoxArray cba = fine.boxArray(); + cba.coarsen(refratio); + cfine.define(cba, fine.DistributionMap(), ncomp, ng); + cfine.setVal(0.0); + cfine.ParallelCopy(crse, 0, 0, ncomp, IntVect(0), ng, crse_geom.periodicity()); + cmf = & cfine; + } + + bool isEB = fine.hasEBFabFactory(); + ignore_unused(isEB); + +#ifdef AMREX_USE_EB + auto factory = dynamic_cast(&(fine.Factory())); + const FabArray* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr; +#endif + + MFItInfo mfi_info; + if (Gpu::notInLaunchRegion()) mfi_info.EnableTiling().SetDynamic(true); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(fine, mfi_info); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const auto& ff = fine.array(mfi); + const auto& cc = cmf->array(mfi); +#ifdef AMREX_USE_EB + bool call_lincc; + if (isEB) + { + const auto& flag = (*flags)[mfi]; + if (flag.getType(amrex::grow(bx,1)) == FabType::regular) { + call_lincc = true; + } else { + Array4 const& flg = flag.const_array(); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, + { + mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp); + }); + + call_lincc = false; + } + } + else + { + call_lincc = true; + } +#else + const bool call_lincc = true; +#endif + if (call_lincc) + { +#if (AMREX_SPACEDIM == 3) + if (hasHiddenDimension()) { + Box const& bx_2d = compactify(bx); + auto const& ff_2d = compactify(ff); + auto const& cc_2d = compactify(cc); + AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx_2d, tbx, + { + TwoD::mlmg_lin_cc_interp_r2(tbx, ff_2d, cc_2d, ncomp); + }); + } else +#endif + { + AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, + { + mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp); + }); + } + } + } +} + +void +MLCellLinOp::AnyAverageDownAndSync (Vector& sol) const +{ + AMREX_ASSERT(sol[0].is()); + + int ncomp = getNComp(); + for (int falev = NAMRLevels()-1; falev > 0; --falev) + { +#ifdef AMREX_USE_EB + amrex::EB_average_down(sol[falev ].get(), + sol[falev-1].get(), 0, ncomp, AMRRefRatio(falev-1)); +#else + amrex::average_down(sol[falev ].get(), + sol[falev-1].get(), 0, ncomp, AMRRefRatio(falev-1)); +#endif + } +} + +void +MLCellLinOp::fixSolvabilityByOffset (int amrlev, int mglev, Any& a_rhs, + Vector const& offset) const +{ + amrex::ignore_unused(amrlev, mglev); + AMREX_ASSERT(a_rhs.is()); + auto& rhs = a_rhs.get(); + + const int ncomp = getNComp(); + for (int c = 0; c < ncomp; ++c) { + rhs.plus(-offset[c], c, 1); + } +#ifdef AMREX_USE_EB + if (rhs.hasEBFabFactory()) { + Vector val(ncomp, 0.0_rt); + amrex::EB_set_covered(rhs, 0, ncomp, val); + } +#endif +} } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H index 1215eda1f6c..1fdf2c60146 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.H @@ -72,7 +72,7 @@ public: virtual std::unique_ptr > makeFactory (int amrlev, int mglev) const final override; - virtual void scaleRHS (int amrlev, MultiFab& rhs) const final; + virtual void scaleRHS (int amrlev, Any& rhs) const final; #endif diff --git a/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp b/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp index cfa7595b515..62a7c3af282 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLEBNodeFDLaplacian.cpp @@ -318,8 +318,11 @@ MLEBNodeFDLaplacian::prepareForSolve () #ifdef AMREX_USE_EB void -MLEBNodeFDLaplacian::scaleRHS (int amrlev, MultiFab& rhs) const +MLEBNodeFDLaplacian::scaleRHS (int amrlev, Any& a_rhs) const { + AMREX_ASSERT(a_rhs.is()); + auto& rhs = a_rhs.get(); + auto const& dmask = *m_dirichlet_mask[amrlev][0]; auto factory = dynamic_cast(m_factory[amrlev][0].get()); auto const& edgecent = factory->getEdgeCent(); @@ -634,19 +637,19 @@ MLEBNodeFDLaplacian::compGrad (int amrlev, const Array #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) void -MLEBNodeFDLaplacian::fillIJMatrix (MFIter const& mfi, - Array4 const& gid, - Array4 const& lid, - HypreNodeLap::Int* const ncols, - HypreNodeLap::Int* const cols, - Real* const mat) const +MLEBNodeFDLaplacian::fillIJMatrix (MFIter const& /*mfi*/, + Array4 const& /*gid*/, + Array4 const& /*lid*/, + HypreNodeLap::Int* const /*ncols*/, + HypreNodeLap::Int* const /*cols*/, + Real* const /*mat*/) const { amrex::Abort("MLEBNodeFDLaplacian::fillIJMatrix: todo"); } void -MLEBNodeFDLaplacian::fillRHS (MFIter const& mfi, Array4 const& lid, - Real* const rhs, Array4 const& bfab) const +MLEBNodeFDLaplacian::fillRHS (MFIter const& /*mfi*/, Array4 const& /*lid*/, + Real* const /*rhs*/, Array4 const& /*bfab*/) const { amrex::Abort("MLEBNodeFDLaplacian::fillRHS: todo"); } diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H index f744c96e059..b42df0e4a42 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H @@ -2,6 +2,7 @@ #define AMREX_ML_LINOP_H_ #include +#include #include #include #include @@ -199,6 +200,7 @@ public: virtual void restriction (int amrlev, int cmglev, MultiFab& crse, MultiFab& fine) const = 0; virtual void interpolation (int amrlev, int fmglev, MultiFab& fine, const MultiFab& crse) const = 0; + virtual void interpAssign (int amrlev, int fmglev, MultiFab& fine, MultiFab& crse) const = 0; virtual void averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, const MultiFab& fine_sol, const MultiFab& fine_rhs) = 0; @@ -223,25 +225,24 @@ public: virtual void compGrad (int amrlev, const Array& grad, MultiFab& sol, Location loc) const = 0; - virtual void applyMetricTerm (int amrlev, int mglev, MultiFab& rhs) const = 0; + virtual void applyMetricTerm (int amrlev, int mglev, Any& rhs) const = 0; virtual void unapplyMetricTerm (int amrlev, int mglev, MultiFab& rhs) const = 0; virtual void fillSolutionBC (int amrlev, MultiFab& sol, const MultiFab* crse_bcdata=nullptr) = 0; - virtual void unimposeNeumannBC (int /*amrlev*/, MultiFab& /*rhs*/) const {} // only nodal solver might need it - virtual void applyInhomogNeumannTerm (int /*amrlev*/, MultiFab& /*rhs*/) const {} - virtual void applyOverset (int /*amlev*/, MultiFab& /*rhs*/) const {} - virtual void scaleRHS (int /*amrlev*/, MultiFab& /*rhs*/) const {} - virtual Real getSolvabilityOffset (int /*amrlev*/, int /*mglev*/, MultiFab const& /*rhs*/) const { return 0._rt; } // Only nodal solvers need it - virtual void fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MultiFab& /*rhs*/, Real /*offset*/) const {} // Only nodal solvers need it + virtual void unimposeNeumannBC (int /*amrlev*/, Any& /*rhs*/) const {} // only nodal solver might need it + virtual void applyInhomogNeumannTerm (int /*amrlev*/, Any& /*rhs*/) const {} + virtual void applyOverset (int /*amlev*/, Any& /*rhs*/) const {} + virtual void scaleRHS (int /*amrlev*/, Any& /*rhs*/) const {} + virtual Vector getSolvabilityOffset (int /*amrlev*/, int /*mglev*/, + Any const& /*rhs*/) const = 0; + virtual void fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, Any& /*rhs*/, + Vector const& /*offset*/) const = 0; virtual void prepareForSolve () = 0; virtual bool isSingular (int amrlev) const = 0; virtual bool isBottomSingular () const = 0; virtual Real xdoty (int amrlev, int mglev, const MultiFab& x, const MultiFab& y, bool local) const = 0; - virtual void fixUpResidualMask (int /*amrlev*/, iMultiFab& /*resmsk*/) { } - virtual void nodalSync (int /*amrlev*/, int /*mglev*/, MultiFab& /*mf*/) const {} - virtual std::unique_ptr makeNLinOp (int grid_size) const = 0; virtual void getFluxes (const Vector >& /*a_flux*/, @@ -283,6 +284,62 @@ public: virtual void copyNSolveSolution (MultiFab&, MultiFab const&) const {} + virtual Any AnyMake (int amrlev, int mglev, int ncomp, IntVect const& ng) const; + virtual Any AnyMakeCoarseMG (int amrlev, int mglev, int ncomp, IntVect const& ng) const; + virtual Any AnyMakeCoarseAmr (int famrlev, int ncomp, IntVect const& ng) const; + virtual Any AnyMakeAlias (Any const& a) const; + virtual IntVect AnyGrowVect (Any const& a) const; + virtual void AnySetBndry (Any& a, Real v, int scomp, int ncomp) const; + virtual void AnyCopy (Any& dst, Any const& src, int scomp, int dcomp, int ncomp, + IntVect const& ng) const; + virtual void AnyAdd (Any& dst, Any const& src, int scomp, int dcomp, int ncomp, + IntVect const& ng) const; + virtual void AnySetToZero (Any& a) const; +#ifdef AMREX_USE_EB + virtual void AnySetCovered (Any& a, int scomp, int ncomp, Vector const& val) const; +#endif + virtual void AnyParallelCopy (Any& dst, Any const& src, + int src_comp, int dst_comp, int ncomp, + IntVect const& src_nghost, IntVect const& dst_nghost, + Periodicity const& period = Periodicity::NonPeriodic()); + virtual void AnyFillBoundary (Any& a, int scomp, int ncomp, IntVect const& nghost, + Periodicity const& period = Periodicity::NonPeriodic()); + + virtual Real AnyNormInf (Any& a) const; + + virtual Real AnyNormInfMask (int amrlev, Any const& a, bool local) const = 0; + + virtual void AnySolutionResidual (int amrlev, Any& resid, Any& x, Any const& b, + Any const* crse_bcdata = nullptr); + virtual void AnyCorrectionResidual (int amrlev, int mglev, Any& resid, Any& x, + const Any& b, BCMode bc_mode, + const Any* crse_bcdata=nullptr); + virtual void AnyReflux (int crse_amrlev, + Any& res, const Any& crse_sol, const Any& crse_rhs, + Any& fine_res, Any& fine_sol, const Any& fine_rhs); + + virtual void AnyAvgDownResAmr (int clev, Any& cres, Any const& fres) const = 0; + virtual void AnyAvgDownResMG (int clev, Any& cres, Any const& fres) const; + + virtual void AnySmooth (int amrlev, int mglev, Any& sol, const Any& rhs, + bool skip_fillboundary=false) const; + + virtual void AnyRestriction (int amrlev, int cmglev, Any& crse, Any& fine) const; + + virtual void AnyInterpolationMG (int amrlev, int fmglev, Any& fine, const Any& crse) const; + virtual void AnyInterpAssignMG (int amrlev, int fmglev, Any& fine, Any& crse) const; + virtual void AnyInterpolationAmr (int famrlev, Any& fine, const Any& crse, + IntVect const& /*nghost*/) const = 0; + + virtual void AnyAverageDownSolutionRHS (int camrlev, Any& crse_sol, Any& crse_rhs, + const Any& fine_sol, const Any& fine_rhs); + + virtual void AnyAverageDownAndSync (Vector& sol) const = 0; + + Real MFNormInf (MultiFab const& mf, iMultiFab const* fine_mask, bool local) const; + + bool isMFIterSafe (int amrlev, int mglev1, int mglev2) const; + protected: static constexpr int mg_coarsen_ratio = 2; @@ -401,7 +458,7 @@ protected: bool isCellCentered () const noexcept { return m_ixtype == 0; } - virtual void make (Vector >& mf, int nc, IntVect const& ng) const; + void make (Vector >& mf, int nc, IntVect const& ng) const; virtual std::unique_ptr > makeFactory (int /*amrlev*/, int /*mglev*/) const { return std::make_unique(); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.cpp index 9c6ccc8ce05..b33a7ef9a0d 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.cpp @@ -4,10 +4,12 @@ #include #include #include +#include #ifdef AMREX_USE_EB #include #include +#include #endif #ifdef AMREX_USE_PETSC @@ -544,7 +546,7 @@ MLLinOp::defineBC () } void -MLLinOp::make (Vector >& mf, int nc, IntVect const& ng) const +MLLinOp::make (Vector >& mf, int nc, IntVect const& ng) const { mf.clear(); mf.resize(m_num_amr_levels); @@ -553,8 +555,7 @@ MLLinOp::make (Vector >& mf, int nc, IntVect const& ng) const mf[alev].resize(m_num_mg_levels[alev]); for (int mlev = 0; mlev < m_num_mg_levels[alev]; ++mlev) { - const auto& ba = amrex::convert(m_grids[alev][mlev], m_ixtype); - mf[alev][mlev].define(ba, m_dmap[alev][mlev], nc, ng, MFInfo(), *m_factory[alev][mlev]); + mf[alev][mlev] = AnyMake(alev, mlev, nc, ng); } } } @@ -895,6 +896,289 @@ MLLinOp::resizeMultiGrid (int new_size) } } +Any +MLLinOp::AnyMake (int amrlev, int mglev, int ncomp, IntVect const& ng) const +{ + return Any(MultiFab(amrex::convert(m_grids[amrlev][mglev], m_ixtype), + m_dmap[amrlev][mglev], ncomp, ng, MFInfo(), + *m_factory[amrlev][mglev])); +} + +Any +MLLinOp::AnyMakeCoarseMG (int amrlev, int mglev, int ncomp, 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 Any(MultiFab(cba, m_dmap[amrlev][mglev], ncomp, ng)); +} + +Any +MLLinOp::AnyMakeCoarseAmr (int famrlev, int ncomp, IntVect const& ng) const +{ + BoxArray cba = m_grids[famrlev][0]; + IntVect ratio(AMRRefRatio(famrlev-1)); + cba.coarsen(ratio); + cba.convert(m_ixtype); + return Any(MultiFab(cba, m_dmap[famrlev][0], ncomp, ng)); +} + +Any +MLLinOp::AnyMakeAlias (Any const& a) const +{ + AMREX_ASSERT(a.is()); + MultiFab const& mf = a.get(); + return Any(MultiFab(mf, amrex::make_alias, 0, mf.nComp())); +} + +IntVect +MLLinOp::AnyGrowVect (Any const& a) const +{ + AMREX_ASSERT(a.is()); + MultiFab const& mf = a.get(); + return mf.nGrowVect(); +} + +void +MLLinOp::AnySetToZero (Any& a) const +{ + AMREX_ASSERT(a.is()); + MultiFab& mf = a.get(); + mf.setVal(0._rt); +} + +void +MLLinOp::AnySetBndry (Any& a, Real v, int scomp, int ncomp) const +{ + AMREX_ASSERT(a.is()); + MultiFab& mf = a.get(); + mf.setBndry(v, scomp, ncomp); +} + +void +MLLinOp::AnyCopy (Any& dst, Any const& src, int scomp, int dcomp, int ncomp, + IntVect const& ng) const +{ + AMREX_ASSERT(dst.is() && src.is()); + MultiFab& dmf = dst.get(); + MultiFab const& smf = src.get(); + MultiFab::Copy(dmf, smf, scomp, dcomp, ncomp, ng); +} + +void +MLLinOp::AnyAdd (Any& dst, Any const& src, int scomp, int dcomp, int ncomp, + IntVect const& ng) const +{ + AMREX_ASSERT(dst.is() && src.is()); + MultiFab& dmf = dst.get(); + MultiFab const& smf = src.get(); + MultiFab::Add(dmf, smf, scomp, dcomp, ncomp, ng); +} + +void +MLLinOp::AnyAverageDownSolutionRHS (int camrlev, Any& a_crse_sol, Any& a_crse_rhs, + const Any& a_fine_sol, const Any& a_fine_rhs) +{ + AMREX_ASSERT(a_crse_sol.is() && + a_crse_rhs.is() && + a_fine_sol.is() && + a_fine_rhs.is()); + auto& crse_sol = a_crse_sol.get(); + auto& crse_rhs = a_crse_rhs.get(); + auto& fine_sol = a_fine_sol.get(); + auto& fine_rhs = a_fine_rhs.get(); + averageDownSolutionRHS(camrlev, crse_sol, crse_rhs, fine_sol, fine_rhs); +} + +#ifdef AMREX_USE_EB +void +MLLinOp::AnySetCovered (Any& a, int scomp, int ncomp, Vector const& val) const +{ + AMREX_ASSERT(a.is()); + auto& mf = a.get(); + EB_set_covered(mf, scomp, ncomp, val); +} +#endif + +void +MLLinOp::AnyParallelCopy (Any& dst, Any const& src, + int src_comp, int dst_comp, int ncomp, + IntVect const& src_nghost, IntVect const& dst_nghost, + Periodicity const& period) +{ + AMREX_ASSERT(dst.is()); + MultiFab& dmf = dst.get(); + MultiFab const& smf = src.get(); + dmf.ParallelCopy(smf, src_comp, dst_comp, ncomp, src_nghost, dst_nghost, period); +} + +void +MLLinOp::AnyFillBoundary (Any& a, int scomp, int ncomp, IntVect const& nghost, + Periodicity const& period) +{ + AMREX_ASSERT(a.is()); + MultiFab& mf = a.get(); + mf.FillBoundary(scomp, ncomp, nghost, period); +} + + +Real +MLLinOp::AnyNormInf (Any& a) const +{ + AMREX_ASSERT(a.is()); + return a.get().norminf(); +} + +void +MLLinOp::AnySolutionResidual (int amrlev, Any& resid, Any& x, Any const& b, + Any const* crse_bcdata) +{ + AMREX_ASSERT(x.is()); + solutionResidual(amrlev, resid.get(), x.get(), b.get(), + (crse_bcdata) ? &(crse_bcdata->get()) : nullptr); +} + +void +MLLinOp::AnyCorrectionResidual (int amrlev, int mglev, Any& resid, Any& x, const Any& b, + BCMode bc_mode, const Any* crse_bcdata) +{ + AMREX_ASSERT(x.is()); + correctionResidual(amrlev, mglev, resid.get(), x.get(), + b.get(), bc_mode, + (crse_bcdata) ? &(crse_bcdata->get()) : nullptr); +} + +void +MLLinOp::AnyReflux (int clev, Any& res, const Any& crse_sol, const Any& crse_rhs, + Any& fine_res, Any& fine_sol, const Any& fine_rhs) +{ + AMREX_ASSERT(res.is()); + reflux(clev,res.get(), crse_sol.get(), crse_rhs.get(), + fine_res.get(), fine_sol.get(), fine_rhs.get()); +} + +Real +MLLinOp::MFNormInf (MultiFab const& mf, iMultiFab const* fine_mask, bool local) const +{ + const int ncomp = getNComp(); + Real norm = 0._rt; + + if (fine_mask == nullptr) { +#ifdef AMREX_USE_GPU + if (Gpu::inLaunchRegion()) { + auto const& ma = mf.const_arrays(); + norm = ParReduce(TypeList{}, TypeList{}, + mf, IntVect(0), ncomp, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) + -> GpuTuple + { + return amrex::Math::abs(ma[box_no](i,j,k,n)); + }); + } else +#endif + { +#ifdef AMREX_USE_OMP +#pragma omp parallel reduction(max:norm) +#endif + for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) { + Box const& bx = mfi.tilebox(); + auto const& fab = mf.const_array(mfi); + AMREX_LOOP_4D(bx, ncomp, i, j, k, n, + { + norm = std::max(norm, amrex::Math::abs(fab(i,j,k,n))); + }); + } + } + } else { +#ifdef AMREX_USE_GPU + if (Gpu::inLaunchRegion()) { + auto const& ma = mf.const_arrays(); + auto const& mask_ma = fine_mask->const_arrays(); + norm = ParReduce(TypeList{}, TypeList{}, + mf, IntVect(0), ncomp, + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k, int n) + -> GpuTuple + { + if (mask_ma[box_no](i,j,k)) { + return amrex::Math::abs(ma[box_no](i,j,k,n)); + } else { + return Real(0.0); + } + }); + } else +#endif + { +#ifdef AMREX_USE_OMP +#pragma omp parallel reduction(max:norm) +#endif + for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) { + Box const& bx = mfi.tilebox(); + auto const& fab = mf.const_array(mfi); + auto const& mask = fine_mask->const_array(mfi); + AMREX_LOOP_4D(bx, ncomp, i, j, k, n, + { + if (mask(i,j,k)) { + norm = std::max(norm, amrex::Math::abs(fab(i,j,k,n))); + } + }); + } + } + } + + if (!local) ParallelAllReduce::Max(norm, ParallelContext::CommunicatorSub()); + return norm; +} + +void +MLLinOp::AnyAvgDownResMG (int clev, Any& cres, Any const& fres) const +{ + AMREX_ASSERT(cres.is()); +#ifdef AMREX_USE_EB + amrex::EB_average_down +#else + amrex::average_down +#endif + (fres.get(), cres.get(), 0, getNComp(), + mg_coarsen_ratio_vec[clev-1]); +} + +void +MLLinOp::AnySmooth (int amrlev, int mglev, Any& sol, const Any& rhs, + bool skip_fillboundary) const +{ + AMREX_ASSERT(sol.is() && rhs.is()); + smooth(amrlev, mglev, sol.get(), rhs.get(), skip_fillboundary); +} + +void +MLLinOp::AnyRestriction (int amrlev, int cmglev, Any& crse, Any& fine) const +{ + AMREX_ASSERT(crse.is() && fine.is()); + restriction(amrlev, cmglev, crse.get(), fine.get()); +} + +void +MLLinOp::AnyInterpolationMG (int amrlev, int fmglev, Any& fine, const Any& crse) const +{ + AMREX_ASSERT(crse.is() && fine.is()); + interpolation(amrlev, fmglev, fine.get(), crse.get()); +} + +void +MLLinOp::AnyInterpAssignMG (int amrlev, int fmglev, Any& fine, Any& crse) const +{ + AMREX_ASSERT(crse.is() && fine.is()); + interpAssign(amrlev, fmglev, fine.get(), crse.get()); +} + +bool +MLLinOp::isMFIterSafe (int amrlev, int mglev1, int mglev2) const +{ + return m_dmap[amrlev][mglev1] == m_dmap[amrlev][mglev2] + && BoxArray::SameRefs(m_grids[amrlev][mglev1], m_grids[amrlev][mglev2]); +} + #ifdef AMREX_USE_PETSC std::unique_ptr MLLinOp::makePETSc () const diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H index 32980d74c45..e884f877fbc 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLMG.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.H @@ -36,6 +36,10 @@ public: Real solve (const Vector& a_sol, const Vector& a_rhs, Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr); + // For this version of solve, Any holds MultiFab like objects. + Real solve (Vector& a_sol, const Vector& a_rhs, + Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file = nullptr); + void getGradSolution (const Vector >& a_grad_sol, Location a_loc = Location::FaceCenter); @@ -121,7 +125,7 @@ public: void setHypreStrongThreshold (Real t) noexcept {hypre_strong_threshold = t;} #endif - void prepareForSolve (const Vector& a_sol, const Vector& a_rhs); + void prepareForSolve (Vector& a_sol, const Vector& a_rhs); void prepareForNSolve (); @@ -151,19 +155,16 @@ public: Real MLRhsNormInf (bool local = false); void buildFineMask (); - void averageDownAndSync (); - - void computeVolInv (); void makeSolvable (); - void makeSolvable (int amrlev, int mglev, MultiFab& mf); + void makeSolvable (int amrlev, int mglev, Any& mf); #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) - void bottomSolveWithHypre (MultiFab& x, const MultiFab& b); + void bottomSolveWithHypre (Any& x, const Any& b); #endif - void bottomSolveWithPETSc (MultiFab& x, const MultiFab& b); + void bottomSolveWithPETSc (Any& x, const Any& b); - int bottomSolveWithCG (MultiFab& x, const MultiFab& b, MLCGSolver::Type type); + int bottomSolveWithCG (Any& x, const Any& b, MLCGSolver::Type type); Real getInitRHS () const noexcept { return m_rhsnorm0; } // Initial composite residual @@ -242,26 +243,21 @@ private: * \brief To avoid confusion, terms like sol, cor, rhs, res, ... etc. are * in the frame of the original equation, not the correction form */ - Vector > sol_raii; - Vector sol; //!< alias to argument a_sol - Vector rhs; //!< Copy of original rhs - //! L(sol) = rhs + Vector sol; //!< Might be alias to argument a_sol + Vector rhs; //!< Copy of original rhs + //! L(sol) = rhs + + Vector sol_is_alias; /** * \brief First Vector: Amr levels. 0 is the coarest level * Second Vector: MG levels. 0 is the finest level */ - Vector > res; //! = rhs - L(sol) - Vector > > cor; //!< L(cor) = res - Vector > > cor_hold; - Vector > rescor; //!< = res - L(cor) - //! Residual of the correction form - - Vector > fine_mask; - - Vector > volinv; //!< used by makeSolvable - - Vector > scratch; + Vector > res; //! = rhs - L(sol) + Vector > cor; //!< L(cor) = res + Vector > cor_hold; + Vector > rescor; //!< = res - L(cor) + //! Residual of the correction form enum timer_types { solve_time=0, iter_time, bottom_time, ntimers }; Vector timer; diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.cpp b/Src/LinearSolvers/MLMG/AMReX_MLMG.cpp index 2bdb9222b4b..8a3c84583c3 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLMG.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.cpp @@ -2,7 +2,6 @@ #include #include #include -#include #include #ifdef AMREX_USE_PETSC @@ -51,25 +50,52 @@ MLMG::~MLMG () Real MLMG::solve (const Vector& a_sol, const Vector& a_rhs, Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file) +{ + Vector any_sol(namrlevs); + Vector any_rhs(namrlevs); + for (int lev = 0; lev < namrlevs; ++lev) { + any_sol[lev] = MultiFab(*a_sol[lev], amrex::make_alias, 0, a_sol[lev]->nComp()); + any_rhs[lev] = MultiFab(*a_rhs[lev], amrex::make_alias, 0, a_rhs[lev]->nComp()); + } + return solve(any_sol, any_rhs, a_tol_rel, a_tol_abs, checkpoint_file); +} + +Real +MLMG::solve (Vector& a_sol, const Vector& a_rhs, + Real a_tol_rel, Real a_tol_abs, const char* checkpoint_file) { BL_PROFILE("MLMG::solve()"); if (checkpoint_file != nullptr) { - checkPoint(a_sol, a_rhs, a_tol_rel, a_tol_abs, checkpoint_file); + if (a_sol[0].is()) { + Vector mf_sol(namrlevs); + Vector mf_rhs(namrlevs); + for (int lev = 0; lev < namrlevs; ++lev) { + mf_sol[lev] = &(a_sol[lev].get()); + mf_rhs[lev] = &(a_rhs[lev].get()); + } + checkPoint(mf_sol, mf_rhs, a_tol_rel, a_tol_abs, checkpoint_file); + } else { + amrex::Abort("MLMG::solve: checkpoint not supported for non-MultiFab type"); + } } if (bottom_solver == BottomSolver::Default) { bottom_solver = linop.getDefaultBottomSolver(); } +#if defined(AMREX_USE_HYPRE) || defined(AMREX_USE_PETSC) if (bottom_solver == BottomSolver::hypre || bottom_solver == BottomSolver::petsc) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_sol[0].is(), + "Non-MultiFab type not supported for hypre and petsc"); int mo = linop.getMaxOrder(); - if (a_sol[0]->hasEBFabFactory()) { + if (a_sol[0].get().hasEBFabFactory()) { linop.setMaxOrder(2); } else { linop.setMaxOrder(std::min(3,mo)); // maxorder = 4 not supported } } +#endif bool is_nsolve = linop.m_parent; @@ -200,9 +226,8 @@ MLMG::solve (const Vector& a_sol, const Vector& a_rh } for (int alev = 0; alev < namrlevs; ++alev) { - if (a_sol[alev] != sol[alev]) - { - MultiFab::Copy(*a_sol[alev], *sol[alev], 0, 0, ncomp, ng_back); + if (!sol_is_alias[alev]) { + linop.AnyCopy(a_sol[alev], sol[alev], 0, 0, ncomp, ng_back); } } @@ -230,15 +255,14 @@ void MLMG::oneIter (int iter) BL_PROFILE("MLMG::oneIter()"); int ncomp = linop.getNComp(); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(); for (int alev = finest_amr_lev; alev > 0; --alev) { - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(alev); miniCycle(alev); - MultiFab::Add(*sol[alev], *cor[alev][0], 0, 0, ncomp, nghost); + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(linop.getNGrow(alev)); + linop.AnyAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost); // compute residual for the coarse AMR level computeResWithCrseSolFineCor(alev-1,alev); @@ -250,7 +274,6 @@ void MLMG::oneIter (int iter) // coarsest amr level { - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(0); // enforce solvability if appropriate if (linop.isSingular(0) && linop.getEnforceSingularSolvable()) { @@ -258,24 +281,27 @@ void MLMG::oneIter (int iter) } if (iter < max_fmg_iters) { - mgFcycle (); + mgFcycle(); } else { - mgVcycle (0, 0); + mgVcycle(0, 0); } - MultiFab::Add(*sol[0], *cor[0][0], 0, 0, ncomp, nghost); + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(linop.getNGrow(0)); + linop.AnyAdd(sol[0], cor[0][0], 0, 0, ncomp, nghost); } for (int alev = 1; alev <= finest_amr_lev; ++alev) { - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(alev); // (Fine AMR correction) = I(Coarse AMR correction) interpCorrection(alev); - MultiFab::Add(*sol[alev], *cor[alev][0], 0, 0, ncomp, nghost); + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(linop.getNGrow(alev)); + linop.AnyAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost); if (alev != finest_amr_lev) { - MultiFab::Add(*cor_hold[alev][0], *cor[alev][0], 0, 0, ncomp, nghost); + linop.AnyAdd(cor_hold[alev][0], cor[alev][0], 0, 0, ncomp, nghost); } // Update fine AMR level correction @@ -283,14 +309,14 @@ void MLMG::oneIter (int iter) miniCycle(alev); - MultiFab::Add(*sol[alev], *cor[alev][0], 0, 0, ncomp, nghost); + linop.AnyAdd(sol[alev], cor[alev][0], 0, 0, ncomp, nghost); if (alev != finest_amr_lev) { - MultiFab::Add(*cor[alev][0], *cor_hold[alev][0], 0, 0, ncomp, nghost); + linop.AnyAdd(cor[alev][0], cor_hold[alev][0], 0, 0, ncomp, nghost); } } - averageDownAndSync(); + linop.AnyAverageDownAndSync(sol); } // Compute multi-level Residual (res) up to amrlevmax. @@ -301,11 +327,11 @@ MLMG::computeMLResidual (int amrlevmax) const int mglev = 0; for (int alev = amrlevmax; alev >= 0; --alev) { - const MultiFab* crse_bcdata = (alev > 0) ? sol[alev-1] : nullptr; - linop.solutionResidual(alev, res[alev][mglev], *sol[alev], rhs[alev], crse_bcdata); + const Any* crse_bcdata = (alev > 0) ? &(sol[alev-1]) : nullptr; + linop.AnySolutionResidual(alev, res[alev][mglev], sol[alev], rhs[alev], crse_bcdata); if (alev < finest_amr_lev) { - linop.reflux(alev, res[alev][mglev], *sol[alev], rhs[alev], - res[alev+1][mglev], *sol[alev+1], rhs[alev+1]); + linop.AnyReflux(alev, res[alev][mglev], sol[alev], rhs[alev], + res[alev+1][mglev], sol[alev+1], rhs[alev+1]); } } } @@ -315,16 +341,8 @@ void MLMG::computeResidual (int alev) { BL_PROFILE("MLMG::computeResidual()"); - - MultiFab& x = *sol[alev]; - const MultiFab& b = rhs[alev]; - MultiFab& r = res[alev][0]; - - const MultiFab* crse_bcdata = nullptr; - if (alev > 0) { - crse_bcdata = sol[alev-1]; - } - linop.solutionResidual(alev, r, x, b, crse_bcdata); + const Any* crse_bcdata = (alev > 0) ? &(sol[alev-1]) : nullptr; + linop.AnySolutionResidual(alev, res[alev][0], sol[alev], rhs[alev], crse_bcdata); } // Compute coarse AMR level composite residual with coarse solution and fine correction @@ -334,38 +352,28 @@ MLMG::computeResWithCrseSolFineCor (int calev, int falev) BL_PROFILE("MLMG::computeResWithCrseSolFineCor()"); int ncomp = linop.getNComp(); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) nghost = std::min(linop.getNGrow(falev),linop.getNGrow(calev)); + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(std::min(linop.getNGrow(falev),linop.getNGrow(calev))); - MultiFab& crse_sol = *sol[calev]; - const MultiFab& crse_rhs = rhs[calev]; - MultiFab& crse_res = res[calev][0]; + Any& crse_sol = sol[calev]; + const Any& crse_rhs = rhs[calev]; + Any& crse_res = res[calev][0]; - MultiFab& fine_sol = *sol[falev]; - const MultiFab& fine_rhs = rhs[falev]; - MultiFab& fine_cor = *cor[falev][0]; - MultiFab& fine_res = res[falev][0]; - MultiFab& fine_rescor = rescor[falev][0]; + Any& fine_sol = sol[falev]; + const Any& fine_rhs = rhs[falev]; + Any& fine_cor = cor[falev][0]; + Any& fine_res = res[falev][0]; + Any& fine_rescor = rescor[falev][0]; - const MultiFab* crse_bcdata = nullptr; - if (calev > 0) { - crse_bcdata = sol[calev-1]; - } - linop.solutionResidual(calev, crse_res, crse_sol, crse_rhs, crse_bcdata); + const Any* crse_bcdata = (calev > 0) ? &(sol[calev-1]) : nullptr; + linop.AnySolutionResidual(calev, crse_res, crse_sol, crse_rhs, crse_bcdata); - linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res, BCMode::Homogeneous); - MultiFab::Copy(fine_res, fine_rescor, 0, 0, ncomp, nghost); + linop.AnyCorrectionResidual(falev, 0, fine_rescor, fine_cor, fine_res, BCMode::Homogeneous); + linop.AnyCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost); - linop.reflux(calev, crse_res, crse_sol, crse_rhs, fine_res, fine_sol, fine_rhs); + linop.AnyReflux(calev, crse_res, crse_sol, crse_rhs, fine_res, fine_sol, fine_rhs); - if (linop.isCellCentered()) { - const int amrrr = linop.AMRRefRatio(calev); -#ifdef AMREX_USE_EB - amrex::EB_average_down(fine_res, crse_res, 0, ncomp, amrrr); -#else - amrex::average_down(fine_res, crse_res, 0, ncomp, amrrr); -#endif - } + linop.AnyAvgDownResAmr(calev, crse_res, fine_res); } // Compute fine AMR level residual fine_res = fine_res - L(fine_cor) with coarse providing BC. @@ -375,19 +383,19 @@ MLMG::computeResWithCrseCorFineCor (int falev) BL_PROFILE("MLMG::computeResWithCrseCorFineCor()"); int ncomp = linop.getNComp(); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(falev); + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(linop.getNGrow(falev)); - const MultiFab& crse_cor = *cor[falev-1][0]; + const Any& crse_cor = cor[falev-1][0]; - MultiFab& fine_cor = *cor[falev][0]; - MultiFab& fine_res = res[falev][0]; - MultiFab& fine_rescor = rescor[falev][0]; + Any& fine_cor = cor [falev][0]; + Any& fine_res = res [falev][0]; + Any& fine_rescor = rescor[falev][0]; // fine_rescor = fine_res - L(fine_cor) - linop.correctionResidual(falev, 0, fine_rescor, fine_cor, fine_res, - BCMode::Inhomogeneous, &crse_cor); - MultiFab::Copy(fine_res, fine_rescor, 0, 0, ncomp, nghost); + linop.AnyCorrectionResidual(falev, 0, fine_rescor, fine_cor, fine_res, + BCMode::Inhomogeneous, &crse_cor); + linop.AnyCopy(fine_res, fine_rescor, 0, 0, ncomp, nghost); } void @@ -413,16 +421,16 @@ MLMG::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { - Real norm = res[amrlev][mglev].norm0(); + Real norm = linop.AnyNormInf(res[amrlev][mglev]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " DN: Norm before smooth " << norm << "\n"; } - cor[amrlev][mglev]->setVal(0.0); + linop.AnySetToZero(cor[amrlev][mglev]); bool skip_fillboundary = true; for (int i = 0; i < nu1; ++i) { - linop.smooth(amrlev, mglev, *cor[amrlev][mglev], res[amrlev][mglev], - skip_fillboundary); + linop.AnySmooth(amrlev, mglev, cor[amrlev][mglev], res[amrlev][mglev], + skip_fillboundary); skip_fillboundary = false; } @@ -431,14 +439,13 @@ MLMG::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { - Real norm = rescor[amrlev][mglev].norm0(); + Real norm = linop.AnyNormInf(rescor[amrlev][mglev]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " DN: Norm after smooth " << norm << "\n"; } // res_crse = R(rescor_fine); this provides res/b to the level below - linop.restriction(amrlev, mglev+1, res[amrlev][mglev+1], rescor[amrlev][mglev]); - + linop.AnyRestriction(amrlev, mglev+1, res[amrlev][mglev+1], rescor[amrlev][mglev]); } BL_PROFILE_VAR("MLMG::mgVcycle_bottom", blp_bottom); @@ -446,7 +453,7 @@ MLMG::mgVcycle (int amrlev, int mglev_top) { if (verbose >= 4) { - Real norm = res[amrlev][mglev_bottom].norm0(); + Real norm = linop.AnyNormInf(res[amrlev][mglev_bottom]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " DN: Norm before bottom " << norm << "\n"; } @@ -454,7 +461,7 @@ MLMG::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { computeResOfCorrection(amrlev, mglev_bottom); - Real norm = rescor[amrlev][mglev_bottom].norm0(); + Real norm = linop.AnyNormInf(rescor[amrlev][mglev_bottom]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " UP: Norm after bottom " << norm << "\n"; @@ -464,21 +471,21 @@ MLMG::mgVcycle (int amrlev, int mglev_top) { if (verbose >= 4) { - Real norm = res[amrlev][mglev_bottom].norm0(); + Real norm = linop.AnyNormInf(res[amrlev][mglev_bottom]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " Norm before smooth " << norm << "\n"; } - cor[amrlev][mglev_bottom]->setVal(0.0); + linop.AnySetToZero(cor[amrlev][mglev_bottom]); bool skip_fillboundary = true; for (int i = 0; i < nu1; ++i) { - linop.smooth(amrlev, mglev_bottom, *cor[amrlev][mglev_bottom], res[amrlev][mglev_bottom], - skip_fillboundary); + linop.AnySmooth(amrlev, mglev_bottom, cor[amrlev][mglev_bottom], + res[amrlev][mglev_bottom], skip_fillboundary); skip_fillboundary = false; } if (verbose >= 4) { computeResOfCorrection(amrlev, mglev_bottom); - Real norm = rescor[amrlev][mglev_bottom].norm0(); + Real norm = linop.AnyNormInf(rescor[amrlev][mglev_bottom]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev_bottom << " Norm after smooth " << norm << "\n"; } @@ -493,12 +500,12 @@ MLMG::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { computeResOfCorrection(amrlev, mglev); - Real norm = rescor[amrlev][mglev].norm0(); + Real norm = linop.AnyNormInf(rescor[amrlev][mglev]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " UP: Norm before smooth " << norm << "\n"; } for (int i = 0; i < nu2; ++i) { - linop.smooth(amrlev, mglev, *cor[amrlev][mglev], res[amrlev][mglev]); + linop.AnySmooth(amrlev, mglev, cor[amrlev][mglev], res[amrlev][mglev]); } if (cf_strategy == CFStrategy::ghostnodes) computeResOfCorrection(amrlev, mglev); @@ -506,7 +513,7 @@ MLMG::mgVcycle (int amrlev, int mglev_top) if (verbose >= 4) { computeResOfCorrection(amrlev, mglev); - Real norm = rescor[amrlev][mglev].norm0(); + Real norm = linop.AnyNormInf(rescor[amrlev][mglev]); amrex::Print() << "AT LEVEL " << amrlev << " " << mglev << " UP: Norm after smooth " << norm << "\n"; } @@ -524,18 +531,12 @@ MLMG::mgFcycle () const int amrlev = 0; const int mg_bottom_lev = linop.NMGLevels(amrlev) - 1; const int ncomp = linop.getNComp(); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(amrlev); + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(linop.getNGrow(amrlev)); for (int mglev = 1; mglev <= mg_bottom_lev; ++mglev) { -#ifdef AMREX_USE_EB - amrex::EB_average_down(res[amrlev][mglev-1], res[amrlev][mglev], 0, ncomp, - linop.mg_coarsen_ratio_vec[mglev-1]); -#else - amrex::average_down(res[amrlev][mglev-1], res[amrlev][mglev], 0, ncomp, - linop.mg_coarsen_ratio_vec[mglev-1]); -#endif + linop.AnyAvgDownResMG(mglev, res[amrlev][mglev], res[amrlev][mglev-1]); } bottomSolve(); @@ -543,17 +544,17 @@ MLMG::mgFcycle () for (int mglev = mg_bottom_lev-1; mglev >= 0; --mglev) { // cor_fine = I(cor_crse) - interpCorrection (amrlev, mglev); + interpCorrection(amrlev, mglev); // rescor = res - L(cor) computeResOfCorrection(amrlev, mglev); // res = rescor; this provides b to the vcycle below - MultiFab::Copy(res[amrlev][mglev], rescor[amrlev][mglev], 0,0,ncomp,nghost); + linop.AnyCopy(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); - MultiFab::Add(*cor[amrlev][mglev], *cor_hold[amrlev][mglev], 0, 0, ncomp, nghost); + linop.AnyAdd(cor[amrlev][mglev], cor_hold[amrlev][mglev], 0, 0, ncomp, nghost); } } @@ -564,16 +565,11 @@ MLMG::interpCorrection (int alev) BL_PROFILE("MLMG::interpCorrection_1"); const int ncomp = linop.getNComp(); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(alev); - - const MultiFab& crse_cor = *cor[alev-1][0]; - MultiFab& fine_cor = *cor[alev][0]; + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(linop.getNGrow(alev)); - BoxArray ba = fine_cor.boxArray(); - const int amrrr = linop.AMRRefRatio(alev-1); - IntVect refratio{amrrr}; - ba.coarsen(refratio); + Any const& crse_cor = cor[alev-1][0]; + Any & fine_cor = cor[alev ][0]; const Geometry& crse_geom = linop.Geom(alev-1,0); @@ -584,121 +580,13 @@ MLMG::interpCorrection (int alev) ng_src = linop.getNGrow(alev-1); ng_dst = linop.getNGrow(alev-1); } - MultiFab cfine(ba, fine_cor.DistributionMap(), ncomp, ng_dst); - cfine.setVal(0.0); - cfine.ParallelCopy(crse_cor, 0, 0, ncomp, ng_src, ng_dst, crse_geom.periodicity()); - - bool isEB = fine_cor.hasEBFabFactory(); - ignore_unused(isEB); - -#ifdef AMREX_USE_EB - auto factory = dynamic_cast(&(fine_cor.Factory())); - const FabArray* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr; -#endif - if (linop.isCellCentered()) - { - MFItInfo mfi_info; - if (Gpu::notInLaunchRegion()) mfi_info.EnableTiling().SetDynamic(true); -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(fine_cor, mfi_info); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - Array4 const& ff = fine_cor.array(mfi); - Array4 const& cc = cfine.const_array(mfi); -#ifdef AMREX_USE_EB - bool call_lincc; - if (isEB) - { - const auto& flag = (*flags)[mfi]; - if (flag.getType(amrex::grow(bx,1)) == FabType::regular) { - call_lincc = true; - } else { - Array4 const& flg = flag.const_array(); - switch(refratio[0]) { - case 2: - { - AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, - { - mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp); - }); - break; - } - case 4: - { - AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, - { - mlmg_eb_cc_interp_r<4>(tbx, ff, cc, flg, ncomp); - }); - break; - } - default: - amrex::Abort("mlmg_eb_cc_interp: only refratio 2 and 4 are supported"); - } + Any cfine = linop.AnyMakeCoarseAmr(alev, ncomp, IntVect(ng_dst)); + linop.AnySetToZero(cfine); + linop.AnyParallelCopy(cfine, crse_cor, 0, 0, ncomp, IntVect(ng_src), IntVect(ng_dst), + crse_geom.periodicity()); - call_lincc = false; - } - } - else - { - call_lincc = true; - } -#else - const bool call_lincc = true; -#endif - if (call_lincc) - { - switch(refratio[0]) { - case 2: - { - AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, - { - mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp); - }); - break; - } - case 4: - { - AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, - { - mlmg_lin_cc_interp_r4(tbx, ff, cc, ncomp); - }); - break; - } - default: - amrex::Abort("mlmg_lin_cc_interp: only refratio 2 and 4 are supported"); - } - } - } - } - else - { - AMREX_ALWAYS_ASSERT(amrrr == 2 || amrrr == 4); -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(fine_cor, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box fbx = mfi.tilebox(); - if (cf_strategy == CFStrategy::ghostnodes && nghost >1) fbx.grow(nghost); - Array4 const& ffab = fine_cor.array(mfi); - Array4 const& cfab = cfine.const_array(mfi); - - if (amrrr == 2) { - AMREX_HOST_DEVICE_FOR_4D ( fbx, ncomp, i, j, k, n, - { - mlmg_lin_nd_interp_r2(i,j,k,n,ffab,cfab); - }); - } else { - AMREX_HOST_DEVICE_FOR_4D ( fbx, ncomp, i, j, k, n, - { - mlmg_lin_nd_interp_r4(i,j,k,n,ffab,cfab); - }); - } - } - } + linop.AnyInterpolationAmr(alev, fine_cor, cfine, nghost); } // Interpolate correction between MG levels @@ -709,119 +597,9 @@ MLMG::interpCorrection (int alev, int mglev) { BL_PROFILE("MLMG::interpCorrection_2"); - MultiFab& crse_cor = *cor[alev][mglev+1]; - MultiFab& fine_cor = *cor[alev][mglev ]; - - const int ncomp = linop.getNComp(); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(alev); - - const Geometry& crse_geom = linop.Geom(alev,mglev+1); - const IntVect refratio = (alev > 0) ? IntVect(2) : linop.mg_coarsen_ratio_vec[mglev]; - - MultiFab cfine; - const MultiFab* cmf; - - if (amrex::isMFIterSafe(crse_cor, fine_cor)) - { - crse_cor.FillBoundary(crse_geom.periodicity()); - cmf = &crse_cor; - } - else - { - BoxArray cba = fine_cor.boxArray(); - cba.coarsen(refratio); - IntVect ng = linop.isCellCentered() ? crse_cor.nGrowVect() : IntVect(0); - if (cf_strategy == CFStrategy::ghostnodes) ng = IntVect(nghost); - cfine.define(cba, fine_cor.DistributionMap(), ncomp, ng); - cfine.setVal(0.0); - cfine.ParallelCopy(crse_cor, 0, 0, ncomp, IntVect(0), ng, crse_geom.periodicity()); - cmf = & cfine; - } - - bool isEB = fine_cor.hasEBFabFactory(); - ignore_unused(isEB); - -#ifdef AMREX_USE_EB - auto factory = dynamic_cast(&(fine_cor.Factory())); - const FabArray* flags = (factory) ? &(factory->getMultiEBCellFlagFab()) : nullptr; -#endif - - if (linop.isCellCentered()) - { - MFItInfo mfi_info; - if (Gpu::notInLaunchRegion()) mfi_info.EnableTiling().SetDynamic(true); -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(fine_cor, mfi_info); mfi.isValid(); ++mfi) - { - const Box& bx = mfi.tilebox(); - const auto& ff = fine_cor.array(mfi); - const auto& cc = cmf->array(mfi); -#ifdef AMREX_USE_EB - bool call_lincc; - if (isEB) - { - const auto& flag = (*flags)[mfi]; - if (flag.getType(amrex::grow(bx,1)) == FabType::regular) { - call_lincc = true; - } else { - Array4 const& flg = flag.const_array(); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, - { - mlmg_eb_cc_interp_r<2>(tbx, ff, cc, flg, ncomp); - }); - - call_lincc = false; - } - } - else - { - call_lincc = true; - } -#else - const bool call_lincc = true; -#endif - if (call_lincc) - { -#if (AMREX_SPACEDIM == 3) - if (linop.hasHiddenDimension()) { - Box const& bx_2d = linop.compactify(bx); - auto const& ff_2d = linop.compactify(ff); - auto const& cc_2d = linop.compactify(cc); - AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx_2d, tbx, - { - TwoD::mlmg_lin_cc_interp_r2(tbx, ff_2d, cc_2d, ncomp); - }); - } else -#endif - { - AMREX_LAUNCH_HOST_DEVICE_LAMBDA (bx, tbx, - { - mlmg_lin_cc_interp_r2(tbx, ff, cc, ncomp); - }); - } - } - } - } - else - { -#ifdef AMREX_USE_OMP -#pragma omp parallel if (Gpu::notInLaunchRegion()) -#endif - for (MFIter mfi(fine_cor, TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - const Box& fbx = mfi.tilebox(); - Array4 const& ffab = fine_cor.array(mfi); - Array4 const& cfab = cmf->const_array(mfi); - - AMREX_HOST_DEVICE_FOR_4D ( fbx, ncomp, i, j, k, n, - { - mlmg_lin_nd_interp_r2(i,j,k,n,ffab,cfab); - }); - } - } + Any& crse_cor = cor[alev][mglev+1]; + Any& fine_cor = cor[alev][mglev ]; + linop.AnyInterpAssignMG(alev, mglev, fine_cor, crse_cor); } // (Fine MG level correction) += I(Coarse MG level correction) @@ -832,29 +610,24 @@ MLMG::addInterpCorrection (int alev, int mglev) const int ncomp = linop.getNComp(); - const MultiFab& crse_cor = *cor[alev][mglev+1]; - MultiFab& fine_cor = *cor[alev][mglev ]; + const Any& crse_cor = cor[alev][mglev+1]; + Any& fine_cor = cor[alev][mglev ]; - MultiFab cfine; - const MultiFab* cmf; + Any cfine; + const Any* cany; - if (amrex::isMFIterSafe(crse_cor, fine_cor)) + if (linop.isMFIterSafe(alev, mglev, mglev+1)) { - cmf = &crse_cor; + cany = &crse_cor; } else { - BoxArray cba = fine_cor.boxArray(); - IntVect ratio = (alev > 0) ? IntVect(2) : linop.mg_coarsen_ratio_vec[mglev]; - - cba.coarsen(ratio); - const int ng = 0; - cfine.define(cba, fine_cor.DistributionMap(), ncomp, ng); - cfine.ParallelCopy(crse_cor); - cmf = &cfine; + cfine = linop.AnyMakeCoarseMG(alev, mglev, ncomp, IntVect(0)); + linop.AnyParallelCopy(cfine,crse_cor,0,0,ncomp,IntVect(0),IntVect(0)); + cany = &cfine; } - linop.interpolation(alev, mglev, fine_cor, *cmf); + linop.AnyInterpolationMG(alev, mglev, fine_cor, *cany); } // Compute rescor = res - L(cor) @@ -865,10 +638,10 @@ void MLMG::computeResOfCorrection (int amrlev, int mglev) { BL_PROFILE("MLMG:computeResOfCorrection()"); - MultiFab& x = *cor[amrlev][mglev]; - const MultiFab& b = res[amrlev][mglev]; - MultiFab& r = rescor[amrlev][mglev]; - linop.correctionResidual(amrlev, mglev, r, x, b, BCMode::Homogeneous); + Any & x = cor[amrlev][mglev]; + const Any& b = res[amrlev][mglev]; + Any & r = rescor[amrlev][mglev]; + linop.AnyCorrectionResidual(amrlev, mglev, r, x, b, BCMode::Homogeneous); } // At the true bottom of the coarset AMR level. @@ -894,7 +667,7 @@ MLMG::NSolve (MLMG& a_solver, MultiFab& a_sol, MultiFab& a_rhs) a_sol.setVal(0.0); - MultiFab const& res_bottom = res[0].back(); + MultiFab const& res_bottom = res[0].back().get(); if (BoxArray::SameRefs(a_rhs.boxArray(),res_bottom.boxArray()) && DistributionMapping::SameRefs(a_rhs.DistributionMap(),res_bottom.DistributionMap())) { @@ -906,7 +679,7 @@ MLMG::NSolve (MLMG& a_solver, MultiFab& a_sol, MultiFab& a_rhs) a_solver.solve({&a_sol}, {&a_rhs}, Real(-1.0), Real(-1.0)); - linop.copyNSolveSolution(*cor[0].back(), a_sol); + linop.copyNSolveSolution(cor[0].back().get(), a_sol); } void @@ -924,28 +697,28 @@ MLMG::actualBottomSolve () const int amrlev = 0; const int mglev = linop.NMGLevels(amrlev) - 1; - MultiFab& x = *cor[amrlev][mglev]; - MultiFab& b = res[amrlev][mglev]; + auto& x = cor[amrlev][mglev]; + auto& b = res[amrlev][mglev]; - x.setVal(0.0); + linop.AnySetToZero(x); if (bottom_solver == BottomSolver::smoother) { bool skip_fillboundary = true; for (int i = 0; i < nuf; ++i) { - linop.smooth(amrlev, mglev, x, b, skip_fillboundary); + linop.AnySmooth(amrlev, mglev, x, b, skip_fillboundary); skip_fillboundary = false; } } else { - MultiFab* bottom_b = &b; - MultiFab raii_b; + Any* bottom_b = &b; + Any raii_b; if (linop.isBottomSingular() && linop.getEnforceSingularSolvable()) { - raii_b.define(b.boxArray(), b.DistributionMap(), ncomp, b.nGrowVect(), - MFInfo(), *linop.Factory(amrlev,mglev)); - MultiFab::Copy(raii_b,b,0,0,ncomp,b.nGrowVect()); + const IntVect ng = linop.AnyGrowVect(b); + raii_b = linop.AnyMake(amrlev, mglev, ncomp, ng); + linop.AnyCopy(raii_b, b, 0, 0, ncomp, ng); bottom_b = &raii_b; makeSolvable(amrlev,mglev,*bottom_b); @@ -973,7 +746,7 @@ MLMG::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(0.0); + linop.AnySetToZero(cor[amrlev][mglev]); if (bottom_solver == BottomSolver::cgbicg || bottom_solver == BottomSolver::bicgcg) { if (bottom_solver == BottomSolver::cgbicg) { @@ -983,7 +756,7 @@ MLMG::actualBottomSolve () } ret = bottomSolveWithCG(x, *bottom_b, cg_type); if (ret != 0) { - cor[amrlev][mglev]->setVal(0.0); + linop.AnySetToZero(cor[amrlev][mglev]); } else { // switch permanently if (cg_type == MLCGSolver::Type::CG) { bottom_solver = BottomSolver::cg; @@ -995,7 +768,7 @@ MLMG::actualBottomSolve () } const int n = (ret==0) ? nub : nuf; for (int i = 0; i < n; ++i) { - linop.smooth(amrlev, mglev, x, b); + linop.AnySmooth(amrlev, mglev, x, b); } } } @@ -1006,7 +779,7 @@ MLMG::actualBottomSolve () } int -MLMG::bottomSolveWithCG (MultiFab& x, const MultiFab& b, MLCGSolver::Type type) +MLMG::bottomSolveWithCG (Any& x, const Any& b, MLCGSolver::Type type) { MLCGSolver cg_solver(this, linop); cg_solver.setSolver(type); @@ -1027,37 +800,7 @@ Real MLMG::ResNormInf (int alev, bool local) { BL_PROFILE("MLMG::ResNormInf()"); - const int ncomp = linop.getNComp(); - const int mglev = 0; - Real norm = 0.0; - MultiFab* pmf = &(res[alev][mglev]); -#ifdef AMREX_USE_EB - if (linop.isCellCentered() && scratch[alev]) { - pmf = scratch[alev].get(); - MultiFab::Copy(*pmf, res[alev][mglev], 0, 0, ncomp, 0); - auto factory = dynamic_cast(linop.Factory(alev)); - if (factory) { - const MultiFab& vfrac = factory->getVolFrac(); - for (int n=0; n < ncomp; ++n) { - MultiFab::Multiply(*pmf, vfrac, 0, n, 1, 0); - } - } else { - amrex::Abort("MLMG::ResNormInf: not EB Factory"); - } - } -#endif - for (int n = 0; n < ncomp; n++) - { - Real newnorm = 0.0; - if (fine_mask[alev]) { - newnorm = pmf->norm0(*fine_mask[alev],n,0,true); - } else { - newnorm = pmf->norm0(n,0,true); - } - norm = std::max(norm, newnorm); - } - if (!local) ParallelAllReduce::Max(norm, ParallelContext::CommunicatorSub()); - return norm; + return linop.AnyNormInfMask(alev, res[alev][0], local); } // Computes multi-level masked inf-norm of Residual (res). @@ -1079,66 +822,17 @@ Real MLMG::MLRhsNormInf (bool local) { BL_PROFILE("MLMG::MLRhsNormInf()"); - const int ncomp = linop.getNComp(); - Real r = 0.0; - for (int alev = 0; alev <= finest_amr_lev; ++alev) - { - MultiFab* pmf = &(rhs[alev]); -#ifdef AMREX_USE_EB - if (linop.isCellCentered() && scratch[alev]) { - pmf = scratch[alev].get(); - MultiFab::Copy(*pmf, rhs[alev], 0, 0, ncomp, 0); - auto factory = dynamic_cast(linop.Factory(alev)); - if (factory) { - const MultiFab& vfrac = factory->getVolFrac(); - for (int n=0; n < ncomp; ++n) { - MultiFab::Multiply(*pmf, vfrac, 0, n, 1, 0); - } - } else { - amrex::Abort("MLMG::MLRhsNormInf: not EB Factory"); - } - } -#endif - for (int n=0; nnorm0(*fine_mask[alev],n,0,true)); - } else { - r = std::max(r, pmf->norm0(n,0,true)); - } - } + Real r = 0.0_rt; + for (int alev = 0; alev <= finest_amr_lev; ++alev) { + auto t = linop.AnyNormInfMask(alev, rhs[alev], true); + r = std::max(r, t); } if (!local) ParallelAllReduce::Max(r, ParallelContext::CommunicatorSub()); return r; } void -MLMG::buildFineMask () -{ - BL_PROFILE("MLMG::buildFineMask()"); - - if (!fine_mask.empty()) return; - - fine_mask.clear(); - fine_mask.resize(namrlevs); - - const auto& amrrr = linop.AMRRefRatio(); - for (int alev = 0; alev < finest_amr_lev; ++alev) - { - fine_mask[alev] = std::make_unique - (makeFineMask(rhs[alev], rhs[alev+1], IntVect(0), IntVect(amrrr[alev]), - Periodicity::NonPeriodic(), 1, 0)); - } - - if (!linop.isCellCentered()) { - for (int alev = 0; alev < finest_amr_lev; ++alev) { - linop.fixUpResidualMask(alev, *fine_mask[alev]); - } - } -} - -void -MLMG::prepareForSolve (const Vector& a_sol, const Vector& a_rhs) +MLMG::prepareForSolve (Vector& a_sol, const Vector& a_rhs) { BL_PROFILE("MLMG::prepareForSolve()"); @@ -1171,29 +865,28 @@ MLMG::prepareForSolve (const Vector& a_sol, const VectornGrowVect() == ng_sol) + else if (linop.AnyGrowVect(a_sol[alev]) == ng_sol) { - sol[alev] = a_sol[alev]; - sol[alev]->setBndry(0.0); + sol[alev] = linop.AnyMakeAlias(a_sol[alev]); + linop.AnySetBndry(sol[alev], 0.0, 0, ncomp); + sol_is_alias[alev] = true; } else { if (!solve_called) { - sol_raii[alev] = std::make_unique(a_sol[alev]->boxArray(), - a_sol[alev]->DistributionMap(), - ncomp, ng_sol, MFInfo(), - *linop.Factory(alev)); + sol[alev] = linop.AnyMake(alev, 0, ncomp, ng_sol); } - MultiFab::Copy(*sol_raii[alev], *a_sol[alev], 0, 0, ncomp, 0); - sol_raii[alev]->setBndry(0.0); - sol[alev] = sol_raii[alev].get(); + linop.AnyCopy(sol[alev], a_sol[alev], 0, 0, ncomp, IntVect(0)); + linop.AnySetBndry(sol[alev], 0.0, 0, ncomp); + sol_is_alias[alev] = false; } } @@ -1202,10 +895,9 @@ MLMG::prepareForSolve (const Vector& a_sol, const VectorboxArray(), a_rhs[alev]->DistributionMap(), ncomp, ng_rhs, - MFInfo(), *linop.Factory(alev)); + rhs[alev] = linop.AnyMake(alev, 0, ncomp, ng_rhs); } - MultiFab::Copy(rhs[alev], *a_rhs[alev], 0, 0, ncomp, ng_rhs); + linop.AnyCopy(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]); @@ -1216,21 +908,21 @@ MLMG::prepareForSolve (const Vector& a_sol, const Vector(linop.Factory(alev)); if (factory) { Vector val(ncomp, 0.0); - amrex::EB_set_covered(rhs[alev], 0, ncomp, val); - amrex::EB_set_covered(*sol[alev], 0, ncomp, val); + linop.AnySetCovered(rhs[alev], 0, ncomp, val); + linop.AnySetCovered(sol[alev], 0, ncomp, val); } #endif } for (int falev = finest_amr_lev; falev > 0; --falev) { - linop.averageDownSolutionRHS(falev-1, *sol[falev-1], rhs[falev-1], *sol[falev], rhs[falev]); + linop.AnyAverageDownSolutionRHS(falev-1, sol[falev-1], rhs[falev-1], + sol[falev], rhs[falev]); } // enforce solvability if appropriate if (linop.isSingular(0) && linop.getEnforceSingularSolvable()) { - computeVolInv(); makeSolvable(); } @@ -1245,8 +937,8 @@ MLMG::prepareForSolve (const Vector& a_sol, const Vector& a_sol, const Vector(res[alev][mglev].boxArray(), - res[alev][mglev].DistributionMap(), - ncomp, _ng, MFInfo(), - *linop.Factory(alev,mglev)); + cor[alev][mglev] = linop.AnyMake(alev, mglev, ncomp, _ng); } - cor[alev][mglev]->setVal(0.0); + linop.AnySetToZero(cor[alev][mglev]); } } @@ -1280,12 +969,9 @@ MLMG::prepareForSolve (const Vector& a_sol, const Vector(cor[alev][mglev]->boxArray(), - cor[alev][mglev]->DistributionMap(), - ncomp, _ng, MFInfo(), - *linop.Factory(alev,mglev)); + cor_hold[alev][mglev] = linop.AnyMake(alev, mglev, ncomp, _ng); } - cor_hold[alev][mglev]->setVal(0.0); + linop.AnySetToZero(cor_hold[alev][mglev]); } } for (int alev = 1; alev < finest_amr_lev; ++alev) @@ -1294,31 +980,9 @@ MLMG::prepareForSolve (const Vector& a_sol, const Vector(cor[alev][0]->boxArray(), - cor[alev][0]->DistributionMap(), - ncomp, _ng, MFInfo(), - *linop.Factory(alev,0)); - } - cor_hold[alev][0]->setVal(0.0); - } - - buildFineMask(); - - if (!solve_called) - { - scratch.resize(namrlevs); -#ifdef AMREX_USE_EB - if (linop.isCellCentered()) { - for (int alev=0; alev < namrlevs; ++alev) { - if (rhs[alev].hasEBFabFactory()) { - scratch[alev] = std::make_unique(rhs[alev].boxArray(), - rhs[alev].DistributionMap(), - ncomp, 0, MFInfo(), - *linop.Factory(alev)); - } - } + cor_hold[alev][0] = linop.AnyMake(alev, 0, ncomp, _ng); } -#endif + linop.AnySetToZero(cor_hold[alev][0]); } if (linop.m_parent) { @@ -1379,7 +1043,7 @@ MLMG::getGradSolution (const Vector >& a_grad_so { BL_PROFILE("MLMG::getGradSolution()"); for (int alev = 0; alev <= finest_amr_lev; ++alev) { - linop.compGrad(alev, a_grad_sol[alev], *sol[alev], a_loc); + linop.compGrad(alev, a_grad_sol[alev], sol[alev].get(), a_loc); } } @@ -1392,7 +1056,11 @@ MLMG::getFluxes (const Vector >& a_flux, } AMREX_ASSERT(sol.size() == a_flux.size()); - getFluxes(a_flux, sol, a_loc); + Vector solmf; + for (auto & s : sol) { + solmf.push_back(&(s.get())); + } + getFluxes(a_flux, solmf, a_loc); } void @@ -1413,7 +1081,11 @@ void MLMG::getFluxes (const Vector & a_flux, Location a_loc) { AMREX_ASSERT(sol.size() == a_flux.size()); - getFluxes(a_flux, sol, a_loc); + Vector solmf; + for (auto & s : sol) { + solmf.push_back(&(s.get())); + } + getFluxes(a_flux, solmf, a_loc); } void @@ -1459,7 +1131,11 @@ MLMG::getEBFluxes (const Vector& a_eb_flux) } AMREX_ASSERT(sol.size() == a_eb_flux.size()); - getEBFluxes(a_eb_flux, sol); + Vector solmf; + for (auto & s : sol) { + solmf.push_back(&(s.get())); + } + getEBFluxes(a_eb_flux, solmf); } void @@ -1486,28 +1162,21 @@ MLMG::compResidual (const Vector& a_res, const Vector& a_s if (linop.hasHiddenDimension()) ng_sol[linop.hiddenDirection()] = 0; sol.resize(namrlevs); - sol_raii.resize(namrlevs); + sol_is_alias.resize(namrlevs,true); for (int alev = 0; alev < namrlevs; ++alev) { - if (cf_strategy == CFStrategy::ghostnodes) - { - sol[alev] = a_sol[alev]; - } - else if (a_sol[alev]->nGrowVect() == ng_sol) + if (cf_strategy == CFStrategy::ghostnodes || a_sol[alev]->nGrowVect() == ng_sol) { - sol[alev] = a_sol[alev]; + sol[alev] = linop.AnyMakeAlias(a_sol[alev]); + sol_is_alias[alev] = true; } else { - if (sol_raii[alev] == nullptr) + if (sol_is_alias[alev]) { - sol_raii[alev] = std::make_unique(a_sol[alev]->boxArray(), - a_sol[alev]->DistributionMap(), - ncomp, ng_sol, MFInfo(), - *linop.Factory(alev)); + sol[alev] = linop.AnyMake(alev, 0, ncomp, ng_sol); } - MultiFab::Copy(*sol_raii[alev], *a_sol[alev], 0, 0, ncomp, 0); - sol[alev] = sol_raii[alev].get(); + MultiFab::Copy(sol[alev].get(), *a_sol[alev], 0, 0, ncomp, 0); } } @@ -1521,22 +1190,23 @@ MLMG::compResidual (const Vector& a_res, const Vector& a_s const auto& amrrr = linop.AMRRefRatio(); for (int alev = finest_amr_lev; alev >= 0; --alev) { - const MultiFab* crse_bcdata = (alev > 0) ? sol[alev-1] : nullptr; + const MultiFab* crse_bcdata = (alev > 0) ? &(sol[alev-1].get()) : nullptr; const MultiFab* prhs = a_rhs[alev]; #if (AMREX_SPACEDIM != 3) int nghost = (cf_strategy == CFStrategy::ghostnodes) ? linop.getNGrow(alev) : 0; - MultiFab rhstmp(prhs->boxArray(), prhs->DistributionMap(), ncomp, nghost, - MFInfo(), *linop.Factory(alev)); + Any rhstmp_a(MultiFab(prhs->boxArray(), prhs->DistributionMap(), ncomp, nghost, + MFInfo(), *linop.Factory(alev))); + MultiFab& rhstmp = rhstmp_a.get(); MultiFab::Copy(rhstmp, *prhs, 0, 0, ncomp, nghost); - linop.applyMetricTerm(alev, 0, rhstmp); - linop.unimposeNeumannBC(alev, rhstmp); - linop.applyInhomogNeumannTerm(alev, rhstmp); + linop.applyMetricTerm(alev, 0, rhstmp_a); + linop.unimposeNeumannBC(alev, rhstmp_a); + linop.applyInhomogNeumannTerm(alev, rhstmp_a); prhs = &rhstmp; #endif - linop.solutionResidual(alev, *a_res[alev], *sol[alev], *prhs, crse_bcdata); + linop.solutionResidual(alev, *a_res[alev], sol[alev].get(), *prhs, crse_bcdata); if (alev < finest_amr_lev) { - linop.reflux(alev, *a_res[alev], *sol[alev], *prhs, - *a_res[alev+1], *sol[alev+1], *a_rhs[alev+1]); + linop.reflux(alev, *a_res[alev], sol[alev].get(), *prhs, + *a_res[alev+1], sol[alev+1].get(), *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]); @@ -1604,7 +1274,8 @@ MLMG::apply (const Vector& out, const Vector& a_in) } for (int alev = 0; alev < namrlevs; ++alev) { - linop.applyInhomogNeumannTerm(alev, rh[alev]); + Any a(MultiFab(rh[alev], amrex::make_alias, 0, rh[alev].nComp())); + linop.applyInhomogNeumannTerm(alev, a); } const auto& amrrr = linop.AMRRefRatio(); @@ -1637,215 +1308,45 @@ MLMG::apply (const Vector& out, const Vector& a_in) } } -void -MLMG::averageDownAndSync () -{ - const auto& amrrr = linop.AMRRefRatio(); - - int ncomp = linop.getNComp(); - int nghost = 0; - if (cf_strategy == CFStrategy::ghostnodes) nghost = linop.getNGrow(); - - if (linop.isCellCentered()) - { - for (int falev = finest_amr_lev; falev > 0; --falev) - { -#ifdef AMREX_USE_EB - amrex::EB_average_down(*sol[falev], *sol[falev-1], 0, ncomp, amrrr[falev-1]); -#else - amrex::average_down(*sol[falev], *sol[falev-1], 0, ncomp, amrrr[falev-1]); -#endif - } - } - else - { - linop.nodalSync(finest_amr_lev, 0, *sol[finest_amr_lev]); - - for (int falev = finest_amr_lev; falev > 0; --falev) - { - const auto& fmf = *sol[falev]; - auto& cmf = *sol[falev-1]; - - MultiFab tmpmf(amrex::coarsen(fmf.boxArray(), amrrr[falev-1]), fmf.DistributionMap(), ncomp, nghost); - amrex::average_down(fmf, tmpmf, 0, ncomp, amrrr[falev-1]); - cmf.ParallelCopy(tmpmf, 0, 0, ncomp); - linop.nodalSync(falev-1, 0, cmf); - } - } -} - -void -MLMG::computeVolInv () -{ - if (solve_called) return; - - if (linop.isCellCentered()) - { - volinv.resize(namrlevs); - for (int amrlev = 0; amrlev < namrlevs; ++amrlev) { - volinv[amrlev].resize(linop.NMGLevels(amrlev)); - } - - // We don't need to compute for every level - - auto f = [&] (int amrlev, int mglev) { -#ifdef AMREX_USE_EB - auto factory = dynamic_cast(linop.Factory(amrlev,mglev)); - if (factory) - { - const MultiFab& vfrac = factory->getVolFrac(); - volinv[amrlev][mglev] = vfrac.sum(0,true); - } - else -#endif - { - volinv[amrlev][mglev] - = Real(1.0 / linop.compactify(linop.Geom(amrlev,mglev).Domain()).d_numPts()); - } - }; - - // amrlev = 0, mglev = 0 - f(0,0); - - int mgbottom = linop.NMGLevels(0)-1; - f(0,mgbottom); - -#ifdef AMREX_USE_EB - Real temp1, temp2; - if (rhs[0].hasEBFabFactory()) - { - ParallelAllReduce::Sum({volinv[0][0], volinv[0][mgbottom]}, - ParallelContext::CommunicatorSub()); - temp1 = Real(1.0)/volinv[0][0]; - temp2 = Real(1.0)/volinv[0][mgbottom]; - } - else - { - temp1 = volinv[0][0]; - temp2 = volinv[0][mgbottom]; - } - volinv[0][0] = temp1; - volinv[0][mgbottom] = temp2; -#endif - } -} - void MLMG::makeSolvable () { - const int ncomp = linop.getNComp(); - - if (linop.isCellCentered()) - { - Vector offset(ncomp); -#ifdef AMREX_USE_EB - auto factory = dynamic_cast(linop.Factory(0)); - if (factory) - { - const MultiFab& vfrac = factory->getVolFrac(); - for (int c = 0; c < ncomp; ++c) { - offset[c] = MultiFab::Dot(rhs[0], c, vfrac, 0, 1, 0, true) * volinv[0][0]; - } - } - else -#endif - { - for (int c = 0; c < ncomp; ++c) { - offset[c] = rhs[0].sum(c,true) * volinv[0][0]; - } - } - ParallelAllReduce::Sum(offset.data(), ncomp, ParallelContext::CommunicatorSub()); - if (verbose >= 4) { - for (int c = 0; c < ncomp; ++c) { - amrex::Print() << "MLMG: Subtracting " << offset[c] - << " from rhs component " << c << "\n"; - } - } - for (int alev = 0; alev < namrlevs; ++alev) { - for (int c = 0; c < ncomp; ++c) { - rhs[alev].plus(-offset[c], c, 1); - } -#ifdef AMREX_USE_EB - if (rhs[alev].hasEBFabFactory()) { - Vector val(ncomp, 0.0); - amrex::EB_set_covered(rhs[alev], 0, ncomp, val); - } -#endif + auto const& offset = linop.getSolvabilityOffset(0, 0, rhs[0]); + if (verbose >= 4) { + const int ncomp = linop.getNComp(); + for (int c = 0; c < ncomp; ++c) { + amrex::Print() << "MLMG: Subtracting " << offset[c] << " from rhs component " + << c << "\n"; } } - else - { - AMREX_ASSERT_WITH_MESSAGE(ncomp==1, "ncomp > 1 not supported for singular nodal problem"); - Real offset = linop.getSolvabilityOffset(0, 0, rhs[0]); - if (verbose >= 4) { - amrex::Print() << "MLMG: Subtracting " << offset << " from rhs\n"; - } - for (int alev = 0; alev < namrlevs; ++alev) { - linop.fixSolvabilityByOffset(alev, 0, rhs[alev], offset); - } + for (int alev = 0; alev < namrlevs; ++alev) { + linop.fixSolvabilityByOffset(alev, 0, rhs[alev], offset); } } void -MLMG::makeSolvable (int amrlev, int mglev, MultiFab& mf) +MLMG::makeSolvable (int amrlev, int mglev, Any& mf) { - const int ncomp = linop.getNComp(); - - if (linop.isCellCentered()) - { - Vector offset(ncomp); -#ifdef AMREX_USE_EB - auto factory = dynamic_cast(linop.Factory(amrlev,mglev)); - if (factory) - { - const MultiFab& vfrac = factory->getVolFrac(); - for (int c = 0; c < ncomp; ++c) { - offset[c] = MultiFab::Dot(mf, c, vfrac, 0, 1, 0, true) * volinv[amrlev][mglev]; - } - } - else -#endif - { - for (int c = 0; c < ncomp; ++c) { - offset[c] = mf.sum(c,true) * volinv[amrlev][mglev]; - } - } - - ParallelAllReduce::Sum(offset.data(), ncomp, ParallelContext::CommunicatorSub()); - - if (verbose >= 4) { - for (int c = 0; c < ncomp; ++c) { - amrex::Print() << "MLMG: Subtracting " << offset[c] - << " from mf component c = " << c << "\n"; - } - } - + auto const& offset = linop.getSolvabilityOffset(amrlev, mglev, mf); + if (verbose >= 4) { + const int ncomp = linop.getNComp(); for (int c = 0; c < ncomp; ++c) { - mf.plus(-offset[c], c, 1); - } -#ifdef AMREX_USE_EB - if (mf.hasEBFabFactory()) { - Vector val(ncomp, 0.0); - amrex::EB_set_covered(mf, 0, ncomp, val); - } -#endif - } - else - { - AMREX_ASSERT_WITH_MESSAGE(ncomp==1, "ncomp > 1 not supported for singular nodal problem"); - Real offset = linop.getSolvabilityOffset(amrlev, mglev, mf); - if (verbose >= 4) { - amrex::Print() << "MLMG: Subtracting " << offset << " on level (" << amrlev << ", " - << mglev << ")\n"; + amrex::Print() << "MLMG: Subtracting " << offset[c] + << " from mf component c = " << c + << " on level (" << amrlev << ", " << mglev << ")\n"; } - linop.fixSolvabilityByOffset(amrlev, mglev, mf, offset); } + linop.fixSolvabilityByOffset(amrlev, mglev, mf, offset); } #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) void -MLMG::bottomSolveWithHypre (MultiFab& x, const MultiFab& b) +MLMG::bottomSolveWithHypre (Any& a_x, const Any& a_b) { + AMREX_ASSERT(a_x.is()); + MultiFab& x = a_x.get(); + MultiFab const& b = a_b.get(); + const int amrlev = 0; const int mglev = linop.NMGLevels(amrlev) - 1; @@ -1905,18 +1406,21 @@ MLMG::bottomSolveWithHypre (MultiFab& x, const MultiFab& b) // For precision reasons we enforce that the average of the correction from hypre is 0 if (linop.isSingular(amrlev) && linop.getEnforceSingularSolvable()) { - makeSolvable(amrlev, mglev, x); + makeSolvable(amrlev, mglev, a_x); } } #endif void -MLMG::bottomSolveWithPETSc (MultiFab& x, const MultiFab& b) +MLMG::bottomSolveWithPETSc (Any& a_x, const Any& a_b) { #if !defined(AMREX_USE_PETSC) - amrex::ignore_unused(x,b); + amrex::ignore_unused(a_x,a_b); amrex::Abort("bottomSolveWithPETSc is called without building with PETSc"); #else + AMREX_ASSERT(a_x.is()); + MultiFab& x = a_x.get(); + MultiFab const& b = a_b.get(); const int ncomp = linop.getNComp(); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ncomp == 1, "bottomSolveWithPETSc doesn't work with ncomp > 1"); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H index affe4c73eaf..50f20e22915 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H @@ -116,9 +116,11 @@ public : } virtual void getFluxes (const Vector& a_flux, const Vector& a_sol) const final override; - virtual void unimposeNeumannBC (int amrlev, MultiFab& rhs) const final override; - virtual Real getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const override; - virtual void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Real offset) const override; + virtual void unimposeNeumannBC (int amrlev, Any& rhs) const final override; + virtual Vector getSolvabilityOffset (int amrlev, int mglev, + Any const& rhs) const override; + virtual void fixSolvabilityByOffset (int amrlev, int mglev, Any& rhs, + Vector const& offset) const override; virtual void compGrad (int /*amrlev*/, const Array& /*grad*/, MultiFab& /*sol*/, Location /*loc*/) const final override { diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp index 79358b58898..c0efaed25d6 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp @@ -150,13 +150,16 @@ MLNodeLaplacian::resizeMultiGrid (int new_size) } void -MLNodeLaplacian::unimposeNeumannBC (int amrlev, MultiFab& rhs) const +MLNodeLaplacian::unimposeNeumannBC (int amrlev, Any& a_rhs) const { if (m_coarsening_strategy == CoarseningStrategy::RAP) { const Box& nddom = amrex::surroundingNodes(Geom(amrlev).Domain()); const auto lobc = LoBC(); const auto hibc = HiBC(); + AMREX_ASSERT(a_rhs.is()); + MultiFab& rhs = a_rhs.get(); + MFItInfo mfi_info; if (Gpu::notInLaunchRegion()) mfi_info.EnableTiling().SetDynamic(true); #ifdef AMREX_USE_OMP @@ -171,14 +174,17 @@ MLNodeLaplacian::unimposeNeumannBC (int amrlev, MultiFab& rhs) const } } -Real -MLNodeLaplacian::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const +Vector +MLNodeLaplacian::getSolvabilityOffset (int amrlev, int mglev, Any const& a_rhs) const { amrex::ignore_unused(amrlev); - AMREX_ASSERT(amrlev==0); - AMREX_ASSERT(mglev+1==m_num_mg_levels[0] || mglev==0); + AMREX_ASSERT(amrlev==0 && (mglev+1==m_num_mg_levels[0] || mglev==0)); + AMREX_ASSERT(getNComp() == 1); if (m_coarsening_strategy == CoarseningStrategy::RAP) { + AMREX_ASSERT(a_rhs.is()); + auto const& rhs = a_rhs.get(); + #ifdef AMREX_USE_EB auto factory = dynamic_cast(m_factory[amrlev][0].get()); if (mglev == 0 && factory && !factory->isAllRegular()) { @@ -229,7 +235,7 @@ MLNodeLaplacian::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rh Real s1 = amrex::get<0>(r); Real s2 = amrex::get<1>(r); ParallelAllReduce::Sum({s1,s2}, ParallelContext::CommunicatorSub()); - return s1/s2; + return {s1/s2}; } else #endif { @@ -279,16 +285,21 @@ MLNodeLaplacian::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rh Real s1 = amrex::get<0>(r); Real s2 = amrex::get<1>(r); ParallelAllReduce::Sum({s1,s2}, ParallelContext::CommunicatorSub()); - return s1/s2; + return {s1/s2}; } } else { - return MLNodeLinOp::getSolvabilityOffset(amrlev, mglev, rhs); + return MLNodeLinOp::getSolvabilityOffset(amrlev, mglev, a_rhs); } } void -MLNodeLaplacian::fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Real offset) const +MLNodeLaplacian::fixSolvabilityByOffset (int amrlev, int mglev, Any& a_rhs, + Vector const& a_offset) const { + AMREX_ASSERT(a_rhs.is()); + auto& rhs = a_rhs.get(); + Real offset = a_offset[0]; + if (m_coarsening_strategy == CoarseningStrategy::RAP) { #ifdef AMREX_USE_EB auto factory = dynamic_cast(m_factory[amrlev][0].get()); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H index c46f4a250f2..1b8f9e7bb56 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H @@ -59,7 +59,7 @@ public: amrex::Abort("AMReX_MLNodeLinOp::compGrad::How did we get here?"); } - virtual void applyMetricTerm (int /*amrlev*/, int /*mglev*/, MultiFab& /*rhs*/) const final override {} + virtual void applyMetricTerm (int /*amrlev*/, int /*mglev*/, Any& /*rhs*/) const final override {} virtual void unapplyMetricTerm (int /*amrlev*/, int /*mglev*/, MultiFab& /*rhs*/) const final override {} virtual void fillSolutionBC (int /*amrlev*/, MultiFab& /*sol*/, @@ -67,12 +67,12 @@ public: amrex::Abort("AMReX_MLNodeLinOp::fillSolutionBC::How did we get here?"); } - virtual void applyInhomogNeumannTerm (int amrlev, MultiFab& rhs) const override; + virtual Vector getSolvabilityOffset (int amrlev, int mglev, + Any const& rhs) const override; + virtual void fixSolvabilityByOffset (int amrlev, int mglev, Any& rhs, + Vector const& offset) const override; - virtual Real getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const override; - virtual void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs, Real offset) const override; - - virtual void prepareForSolve () override {} + virtual void prepareForSolve () override; virtual bool isSingular (int amrlev) const override { return (amrlev == 0) ? m_is_bottom_singular : false; } @@ -86,7 +86,7 @@ public: virtual void Fapply (int amrlev, int mglev, MultiFab& out, const MultiFab& in) const = 0; virtual void Fsmooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rsh) const = 0; - virtual void nodalSync (int amrlev, int mglev, MultiFab& mf) const final override; + void nodalSync (int amrlev, int mglev, MultiFab& mf) const; virtual std::unique_ptr makeNLinOp (int /*grid_size*/) const final override { amrex::Abort("MLNodeLinOp::makeNLinOp: N-Solve not supported"); @@ -102,6 +102,19 @@ public: // omask is either 0 or 1. 1 means the node is an unknown. 0 means it's known. void setOversetMask (int amrlev, const iMultiFab& a_omask); + virtual void fixUpResidualMask (int /*amrlev*/, iMultiFab& /*resmsk*/) { } + + virtual Real AnyNormInfMask (int amrlev, Any const& a, bool local) const override; + + virtual void AnyAvgDownResAmr (int, Any&, Any const&) const final override { } + + virtual void AnyInterpolationAmr (int famrlev, Any& fine, const Any& crse, + IntVect const& nghost) const override; + + virtual void AnyAverageDownAndSync (Vector& sol) const override; + + virtual void interpAssign (int amrlev, int fmglev, MultiFab& fine, MultiFab& crse) const override; + #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) virtual std::unique_ptr makeHypreNodeLap( int bottom_verbose, @@ -139,6 +152,8 @@ protected: MultiFab m_bottom_dot_mask; MultiFab m_coarse_dot_mask; + Vector > m_norm_fine_mask; + #ifdef AMREX_USE_EB CoarseningStrategy m_coarsening_strategy = CoarseningStrategy::RAP; #else diff --git a/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.cpp b/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.cpp index baf0f5edb42..b5173b71f5f 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.cpp +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.cpp @@ -1,6 +1,7 @@ #include #include +#include #include #ifdef AMREX_USE_OMP @@ -83,6 +84,22 @@ MLNodeLinOp::define (const Vector& a_geom, m_has_fine_bndry[amrlev] = std::make_unique >(m_grids[amrlev][0], m_dmap[amrlev][0]); } + + m_norm_fine_mask.resize(m_num_amr_levels-1); + for (int amrlev = 0; amrlev < m_num_amr_levels-1; ++amrlev) { + m_norm_fine_mask[amrlev] = std::make_unique + (makeFineMask(amrex::convert(m_grids[amrlev][0], IntVect(1)), m_dmap[amrlev][0], + amrex::convert(m_grids[amrlev+1][0], IntVect(1)), + IntVect(m_amr_ref_ratio[amrlev]), 1, 0)); + } +} + +void +MLNodeLinOp::prepareForSolve () +{ + for (int amrlev = 0; amrlev < m_num_amr_levels-1; ++amrlev) { + fixUpResidualMask(amrlev, *m_norm_fine_mask[amrlev]); + } } std::unique_ptr @@ -177,17 +194,16 @@ MLNodeLinOp::xdoty (int amrlev, int mglev, const MultiFab& x, const MultiFab& y, return result; } -void -MLNodeLinOp::applyInhomogNeumannTerm (int /*amrlev*/, MultiFab& /*rhs*/) const -{ -} - -Real -MLNodeLinOp::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) const +Vector +MLNodeLinOp::getSolvabilityOffset (int amrlev, int mglev, Any const& a_rhs) const { amrex::ignore_unused(amrlev); - AMREX_ASSERT(amrlev==0); - AMREX_ASSERT(mglev+1==m_num_mg_levels[0] || mglev==0); + AMREX_ASSERT(amrlev==0 && (mglev+1==m_num_mg_levels[0] || mglev==0)); + AMREX_ASSERT(getNComp() == 1); + + AMREX_ASSERT(a_rhs.is()); + auto const& rhs = a_rhs.get(); + const auto& mask = (mglev+1 == m_num_mg_levels[0]) ? m_bottom_dot_mask : m_coarse_dot_mask; const auto& mask_ma = mask.const_arrays(); const auto& rhs_ma = rhs.const_arrays(); @@ -203,13 +219,16 @@ MLNodeLinOp::getSolvabilityOffset (int amrlev, int mglev, MultiFab const& rhs) c Real s1 = amrex::get<0>(r); Real s2 = amrex::get<1>(r); ParallelAllReduce::Sum({s1,s2}, ParallelContext::CommunicatorSub()); - return s1/s2; + return {s1/s2}; } void -MLNodeLinOp::fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, MultiFab& rhs, Real offset) const +MLNodeLinOp::fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, Any& a_rhs, + Vector const& offset) const { - rhs.plus(-offset, 0, 1); + AMREX_ASSERT(a_rhs.is()); + auto& rhs = a_rhs.get(); + rhs.plus(-offset[0], 0, 1); } namespace { @@ -448,6 +467,119 @@ MLNodeLinOp::resizeMultiGrid (int new_size) MLLinOp::resizeMultiGrid(new_size); } +Real +MLNodeLinOp::AnyNormInfMask (int amrlev, Any const& a, bool local) const +{ + AMREX_ASSERT(a.is()); + auto& mf = a.get(); + + const int finest_level = NAMRLevels() - 1; + iMultiFab const* fine_mask = (amrlev == finest_level) + ? nullptr : m_norm_fine_mask[amrlev].get(); + return MFNormInf(mf, fine_mask, local); +} + +void +MLNodeLinOp::AnyInterpolationAmr (int famrlev, Any& a_fine, const Any& a_crse, + IntVect const& nghost) const +{ + AMREX_ASSERT(a_fine.is()); + MultiFab& fine = a_fine.get(); + MultiFab const& crse = a_crse.get(); + + const int ncomp = getNComp(); + const int refratio = AMRRefRatio(famrlev-1); + + AMREX_ALWAYS_ASSERT(refratio == 2 || refratio == 4); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(fine, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box fbx = mfi.tilebox(); + fbx.grow(nghost); + Array4 const& ffab = fine.array(mfi); + Array4 const& cfab = crse.const_array(mfi); + + if (refratio == 2) { + AMREX_HOST_DEVICE_FOR_4D ( fbx, ncomp, i, j, k, n, + { + mlmg_lin_nd_interp_r2(i,j,k,n,ffab,cfab); + }); + } else { + AMREX_HOST_DEVICE_FOR_4D ( fbx, ncomp, i, j, k, n, + { + mlmg_lin_nd_interp_r4(i,j,k,n,ffab,cfab); + }); + } + } +} + +void +MLNodeLinOp::AnyAverageDownAndSync (Vector& sol) const +{ + AMREX_ASSERT(sol[0].is()); + + const int ncomp = getNComp(); + const int finest_amr_lev = NAMRLevels() - 1; + + nodalSync(finest_amr_lev, 0, sol[finest_amr_lev].get()); + + for (int falev = finest_amr_lev; falev > 0; --falev) + { + const auto& fmf = sol[falev ].get(); + auto& cmf = sol[falev-1].get(); + + auto rr = AMRRefRatio(falev-1); + MultiFab tmpmf(amrex::coarsen(fmf.boxArray(), rr), fmf.DistributionMap(), ncomp, 0); + amrex::average_down(fmf, tmpmf, 0, ncomp, rr); + cmf.ParallelCopy(tmpmf, 0, 0, ncomp); + nodalSync(falev-1, 0, cmf); + } +} + +void +MLNodeLinOp::interpAssign (int amrlev, int fmglev, MultiFab& fine, MultiFab& crse) const +{ + const int ncomp = getNComp(); + + const Geometry& crse_geom = Geom(amrlev,fmglev+1); + const IntVect refratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[fmglev]; + AMREX_ALWAYS_ASSERT(refratio == 2); + + MultiFab cfine; + const MultiFab* cmf; + + if (amrex::isMFIterSafe(crse, fine)) + { + crse.FillBoundary(crse_geom.periodicity()); + cmf = &crse; + } + else + { + BoxArray cba = fine.boxArray(); + cba.coarsen(refratio); + cfine.define(cba, fine.DistributionMap(), ncomp, 0); + cfine.ParallelCopy(crse, 0, 0, ncomp, 0, 0, crse_geom.periodicity()); + cmf = & cfine; + } + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(fine, TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const Box& fbx = mfi.tilebox(); + Array4 const& ffab = fine.array(mfi); + Array4 const& cfab = cmf->const_array(mfi); + + AMREX_HOST_DEVICE_FOR_4D ( fbx, ncomp, i, j, k, n, + { + mlmg_lin_nd_interp_r2(i,j,k,n,ffab,cfab); + }); + } +} + #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1) std::unique_ptr MLNodeLinOp::makeHypreNodeLap (int bottom_verbose, const std::string& options_namespace) const