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/Base/AMReX_BaseFab.H b/Src/Base/AMReX_BaseFab.H index 3a9f5eea018..60f93170337 100644 --- a/Src/Base/AMReX_BaseFab.H +++ b/Src/Base/AMReX_BaseFab.H @@ -260,7 +260,7 @@ public: */ void clear () noexcept; - // Release ownership of memory + //! Release ownership of memory std::unique_ptr release () noexcept; //! Returns how many bytes used diff --git a/Src/Base/AMReX_FabArray.H b/Src/Base/AMReX_FabArray.H index 6eef7caa579..8be30fc8763 100644 --- a/Src/Base/AMReX_FabArray.H +++ b/Src/Base/AMReX_FabArray.H @@ -2848,7 +2848,7 @@ FabArray::SumBoundary_nowait (int scomp, int ncomp, IntVect const& src_ngho FabArray* tmp = new FabArray( boxArray(), DistributionMap(), ncomp, src_nghost, MFInfo(), Factory() ); amrex::Copy(*tmp, *this, scomp, 0, ncomp, src_nghost); - this->setVal(0.0, scomp, ncomp, dst_nghost); + this->setVal(typename FAB::value_type(0), scomp, ncomp, dst_nghost); this->ParallelCopy_nowait(*tmp,0,scomp,ncomp,src_nghost,dst_nghost,period,FabArrayBase::ADD); // All local. Operation complete. 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..457f7565df3 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,9 +135,12 @@ 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; @@ -146,6 +152,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 +228,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..e4c9cef953f 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) @@ -530,18 +540,6 @@ MLCellLinOp::solutionResidual (int amrlev, MultiFab& resid, MultiFab& x, const M MultiFab::Xpay(resid, Real(-1.0), b, 0, 0, ncomp, 0); } -void -MLCellLinOp::fillSolutionBC (int amrlev, MultiFab& sol, const MultiFab* crse_bcdata) -{ - BL_PROFILE("MLCellLinOp::fillSolutionBC()"); - if (crse_bcdata != nullptr) { - updateSolBC(amrlev, *crse_bcdata); - } - const int mglev = 0; - applyBC(amrlev, mglev, sol, BCMode::Inhomogeneous, StateMode::Solution, - m_bndry_sol[amrlev].get()); -} - void MLCellLinOp::correctionResidual (int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b, BCMode bc_mode, const MultiFab* crse_bcdata) @@ -1316,7 +1314,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 +1446,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..f7096b93778 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 @@ -177,10 +178,10 @@ public: * inhomogeneous Neumann BC, the value in leveldata is assumed to be * `d./dx`. */ - virtual void setLevelBC (int amrlev, const MultiFab* levelbcdata, - const MultiFab* robinbc_a = nullptr, - const MultiFab* robinbc_b = nullptr, - const MultiFab* robinbc_f = nullptr) = 0; + virtual void setLevelBC (int /*amrlev*/, const MultiFab* /*levelbcdata*/, + const MultiFab* /*robinbc_a*/ = nullptr, + const MultiFab* /*robinbc_b*/ = nullptr, + const MultiFab* /*robinbc_f*/ = nullptr) {} void setVerbose (int v) noexcept { verbose = v; } @@ -197,52 +198,51 @@ public: virtual bool needsUpdate () const { return false; } virtual void update () {} - 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 averageDownSolutionRHS (int camrlev, MultiFab& crse_sol, MultiFab& crse_rhs, - const MultiFab& fine_sol, const MultiFab& fine_rhs) = 0; + virtual void restriction (int /*amrlev*/, int /*cmglev*/, MultiFab& /*crse*/, MultiFab& /*fine*/) const {} + virtual void interpolation (int /*amrlev*/, int /*fmglev*/, MultiFab& /*fine*/, const MultiFab& /*crse*/) const {} + virtual void interpAssign (int /*amrlev*/, int /*fmglev*/, MultiFab& /*fine*/, MultiFab& /*crse*/) const {} + virtual void averageDownSolutionRHS (int /*camrlev*/, MultiFab& /*crse_sol*/, MultiFab& /*crse_rhs*/, + const MultiFab& /*fine_sol*/, const MultiFab& /*fine_rhs*/) {} - virtual void apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc_mode, - StateMode s_mode, const MLMGBndry* bndry=nullptr) const = 0; - virtual void smooth (int amrlev, int mglev, MultiFab& sol, const MultiFab& rhs, - bool skip_fillboundary=false) const = 0; + virtual void apply (int /*amrlev*/, int /*mglev*/, MultiFab& /*out*/, MultiFab& /*in*/, BCMode /*bc_mode*/, + StateMode /*s_mode*/, const MLMGBndry* /*bndry*/=nullptr) const {} + virtual void smooth (int /*amrlev*/, int /*mglev*/, MultiFab& /*sol*/, const MultiFab& /*rhs*/, + bool /*skip_fillboundary*/=false) const {} // Divide mf by the diagonal component of the operator. Used by bicgstab. virtual void normalize (int /*amrlev*/, int /*mglev*/, MultiFab& /*mf*/) const {} - virtual void solutionResidual (int amrlev, MultiFab& resid, MultiFab& x, const MultiFab& b, - const MultiFab* crse_bcdata=nullptr) = 0; - virtual void correctionResidual (int amrlev, int mglev, MultiFab& resid, MultiFab& x, const MultiFab& b, - BCMode bc_mode, const MultiFab* crse_bcdata=nullptr) = 0; - - virtual void reflux (int crse_amrlev, - MultiFab& res, const MultiFab& crse_sol, const MultiFab& crse_rhs, - MultiFab& fine_res, MultiFab& fine_sol, const MultiFab& fine_rhs) const = 0; - virtual void compFlux (int amrlev, const Array& fluxes, - MultiFab& sol, Location loc) const = 0; - 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 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 solutionResidual (int /*amrlev*/, MultiFab& /*resid*/, MultiFab& /*x*/, const MultiFab& /*b*/, + const MultiFab* /*crse_bcdata*/=nullptr) {} + virtual void correctionResidual (int /*amrlev*/, int /*mglev*/, MultiFab& /*resid*/, MultiFab& /*x*/, const MultiFab& /*b*/, + BCMode /*bc_mode*/, const MultiFab* /*crse_bcdata*/=nullptr) {} + + virtual void reflux (int /*crse_amrlev*/, + MultiFab& /*res*/, const MultiFab& /*crse_sol*/, const MultiFab& /*crse_rhs*/, + MultiFab& /*fine_res*/, MultiFab& /*fine_sol*/, const MultiFab& /*fine_rhs*/) const {} + virtual void compFlux (int /*amrlev*/, const Array& /*fluxes*/, + MultiFab& /*sol*/, Location /*loc*/) const {} + virtual void compGrad (int /*amrlev*/, const Array& /*grad*/, + MultiFab& /*sol*/, Location /*loc*/) const {} + + virtual void applyMetricTerm (int /*amrlev*/, int /*mglev*/, Any& /*rhs*/) const {} + virtual void unapplyMetricTerm (int /*amrlev*/, int /*mglev*/, MultiFab& /*rhs*/) const {} + + 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 { return {}; } + virtual void fixSolvabilityByOffset (int /*amrlev*/, int /*mglev*/, Any& /*rhs*/, + Vector const& /*offset*/) const {} 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 bool isSingular (int /*amrlev*/) const { return false; } + virtual bool isBottomSingular () const { return false; } + virtual Real xdoty (int /*amrlev*/, int /*mglev*/, const MultiFab& /*x*/, const MultiFab& /*y*/, bool /*local*/) const { return 0._rt; } - 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 std::unique_ptr makeNLinOp (int /*grid_size*/) const { return {nullptr}; } virtual void getFluxes (const Vector >& /*a_flux*/, const Vector& /*a_sol*/, @@ -283,6 +283,57 @@ public: virtual void copyNSolveSolution (MultiFab&, MultiFab const&) const {} + virtual Any AnyMake (int amrlev, int mglev, IntVect const& ng) const; + virtual Any AnyMakeCoarseMG (int amrlev, int mglev, IntVect const& ng) const; + virtual Any AnyMakeCoarseAmr (int famrlev, IntVect const& ng) const; + virtual Any AnyMakeAlias (Any const& a) const; + virtual IntVect AnyGrowVect (Any const& a) const; + virtual void AnyCopy (Any& dst, Any const& src, IntVect const& ng) const; + virtual void AnyAdd (Any& dst, Any const& src, IntVect const& ng) const; + virtual void AnySetToZero (Any& a) const; + virtual void AnySetBndryToZero (Any& a) const; +#ifdef AMREX_USE_EB + virtual void AnySetCoveredToZero (Any& a) const; +#endif + virtual void AnyParallelCopy (Any& dst, Any const& src, + IntVect const& src_nghost, IntVect const& dst_nghost, + Periodicity const& period = Periodicity::NonPeriodic()) const; + + 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 +452,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, 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..5f71895320d 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, 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, ng); } } } @@ -895,6 +896,276 @@ MLLinOp::resizeMultiGrid (int new_size) } } +Any +MLLinOp::AnyMake (int amrlev, int mglev, IntVect const& ng) const +{ + return Any(MultiFab(amrex::convert(m_grids[amrlev][mglev], m_ixtype), + m_dmap[amrlev][mglev], getNComp(), ng, MFInfo(), + *m_factory[amrlev][mglev])); +} + +Any +MLLinOp::AnyMakeCoarseMG (int amrlev, int mglev, IntVect const& ng) const +{ + BoxArray cba = m_grids[amrlev][mglev]; + IntVect ratio = (amrlev > 0) ? IntVect(2) : mg_coarsen_ratio_vec[mglev]; + cba.coarsen(ratio); + cba.convert(m_ixtype); + return Any(MultiFab(cba, m_dmap[amrlev][mglev], getNComp(), ng)); +} + +Any +MLLinOp::AnyMakeCoarseAmr (int famrlev, IntVect const& ng) const +{ + BoxArray cba = m_grids[famrlev][0]; + IntVect ratio(AMRRefRatio(famrlev-1)); + cba.coarsen(ratio); + cba.convert(m_ixtype); + return Any(MultiFab(cba, m_dmap[famrlev][0], getNComp(), 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::AnySetBndryToZero (Any& a) const +{ + AMREX_ASSERT(a.is()); + MultiFab& mf = a.get(); + mf.setBndry(0._rt, 0, getNComp()); +} + +#ifdef AMREX_USE_EB +void +MLLinOp::AnySetCoveredToZero (Any& a) const +{ + AMREX_ASSERT(a.is()); + auto& mf = a.get(); + EB_set_covered(mf, 0, getNComp(), 0, 0._rt); +} +#endif + +void +MLLinOp::AnyCopy (Any& dst, Any const& src, IntVect const& ng) const +{ + AMREX_ASSERT(dst.is() && src.is()); + MultiFab& dmf = dst.get(); + MultiFab const& smf = src.get(); + MultiFab::Copy(dmf, smf, 0, 0, getNComp(), ng); +} + +void +MLLinOp::AnyAdd (Any& dst, Any const& src, IntVect const& ng) const +{ + AMREX_ASSERT(dst.is() && src.is()); + MultiFab& dmf = dst.get(); + MultiFab const& smf = src.get(); + MultiFab::Add(dmf, smf, 0, 0, getNComp(), 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); +} + +void +MLLinOp::AnyParallelCopy (Any& dst, Any const& src, + IntVect const& src_nghost, IntVect const& dst_nghost, + Periodicity const& period) const +{ + AMREX_ASSERT(dst.is()); + MultiFab& dmf = dst.get(); + MultiFab const& smf = src.get(); + dmf.ParallelCopy(smf, 0, 0, getNComp(), src_nghost, dst_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_MLLinOp_temp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp_temp.H new file mode 100644 index 00000000000..68d7c836ba5 --- /dev/null +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp_temp.H @@ -0,0 +1,486 @@ +#ifndef AMREX_MLLINOP_TEMP_H_ +#define AMREX_MLLINOP_TEMP_H_ + +//! This is a template for writing your own linear operator class for Ax=b. + +#include + +namespace amrex_temp +{ + +class MLLinOpTemp + : public amrex::MLLinOp +{ +public: + + //! In this example, there are 3 edge based MultiFabs. + using Container = amrex::Array; + + MLLinOpTemp () {} + + virtual ~MLLinOpTemp () {} + + MLLinOpTemp (const MLLinOpTemp&) = delete; + MLLinOpTemp (MLLinOpTemp&&) = delete; + MLLinOpTemp& operator= (const MLLinOpTemp&) = delete; + MLLinOpTemp& operator= (MLLinOpTemp&&) = delete; + + MLLinOpTemp (const amrex::Vector& a_geom, + const amrex::Vector& a_grids, + const amrex::Vector& a_dmap, + const amrex::LPInfo& a_info = amrex::LPInfo(), + const amrex::Vector const*>& a_factory = {}) + { + define(a_geom, a_grids, a_dmap, a_info, a_factory); + } + + void define (const amrex::Vector& a_geom, + const amrex::Vector& a_grids, + const amrex::Vector& a_dmap, + const amrex::LPInfo& a_info = amrex::LPInfo(), + const amrex::Vector const*>& a_factory = {}) + { + amrex::MLLinOp::define(a_geom, a_grids, a_dmap, a_info, a_factory); + } + + /** + * \brief Return the default solver at the bottom of MG cycles. By + * default, MLLinOp uses a BiCGStab solver implemented in + * AMReX::MLCGSolver. However, it only supports a single MultiFab. + * Since our data type is different, we use a smoother instead. In the + * future we can try to generalize MLCGSolver. + */ + virtual amrex::BottomSolver getDefaultBottomSolver () const override { + return amrex::BottomSolver::smoother; + } + + /** + * \brief Make data container (e.g., MultiFabs stored in Any) for given level. + * + * \param amrlev AMR level. Note that the lowest level is always 0. + * \param mglev MG level. Note that mglev+1 is one level coarser than mglev. + * \param ng number of ghost cells. + */ + virtual amrex::Any AnyMake (int amrlev, int mglev, amrex::IntVect const& ng) const override + { + auto const& ba = m_grids[amrlev][mglev]; + auto const& dm = m_dmap [amrlev][mglev]; + auto const& fc = *m_factory[amrlev][mglev]; + return amrex::Any(Container{amrex::MultiFab(amrex::convert(ba,amrex::IntVect(0,1,1)), + dm, 1, ng, amrex::MFInfo(), fc), + amrex::MultiFab(amrex::convert(ba,amrex::IntVect(1,0,1)), + dm, 1, ng, amrex::MFInfo(), fc), + amrex::MultiFab(amrex::convert(ba,amrex::IntVect(1,1,0)), + dm, 1, ng, amrex::MFInfo(), fc)}); + } + + /** + * \brief Make data container with coarsened BoxArray and + * DistributionMapping of the give MG level. + * + * \param amrlev AMR level. Note that the lowest level is always 0. + * \param mglev MG level. The coarser level is mglev+1. + * \param ng number of ghost cells. + */ + virtual amrex::Any AnyMakeCoarseMG (int amrlev, int mglev, amrex::IntVect const& ng) const override + { + auto ratio = (amrlev > 0) ? amrex::IntVect(2) : this->mg_coarsen_ratio_vec[mglev]; + auto const& ba = amrex::coarsen(m_grids[amrlev][mglev], ratio); + auto const& dm = m_dmap[amrlev][mglev]; + return amrex::Any(Container{amrex::MultiFab(amrex::convert(ba,amrex::IntVect(0,1,1)), + dm, 1, ng), + amrex::MultiFab(amrex::convert(ba,amrex::IntVect(1,0,1)), + dm, 1, ng), + amrex::MultiFab(amrex::convert(ba,amrex::IntVect(1,1,0)), + dm, 1, ng)}); + } + + /** + * \brief Make data container with coarsened BoxArray and + * DistributionMapping of the given AMR level. + * + * \param famrlev AMR level. The coarser AMR level is famrlev-1. + * \param ng number of ghost cells. + */ + virtual amrex::Any AnyMakeCoarseAmr (int famrlev, amrex::IntVect const& ng) const override + { + amrex::IntVect ratio(this->AMRRefRatio(famrlev-1)); + auto const& ba = amrex::coarsen(m_grids[famrlev][0], ratio); + auto const& dm = m_dmap[famrlev][0]; + return amrex::Any(Container{amrex::MultiFab(amrex::convert(ba,amrex::IntVect(0,1,1)), + dm, 1, ng), + amrex::MultiFab(amrex::convert(ba,amrex::IntVect(1,0,1)), + dm, 1, ng), + amrex::MultiFab(amrex::convert(ba,amrex::IntVect(1,1,0)), + dm, 1, ng)}); + } + + /** + * \brief Make an alias of the given Any without deepcopying. + * + * \param a an Any object. + */ + virtual amrex::Any AnyMakeAlias (amrex::Any const& a) const override + { + auto const& rhs = a.get(); + return amrex::Any(Container{amrex::MultiFab(rhs[0], amrex::make_alias, 0, 1), + amrex::MultiFab(rhs[1], amrex::make_alias, 0, 1), + amrex::MultiFab(rhs[2], amrex::make_alias, 0, 1)}); + } + + /** + * \brief Retuen the number of ghost cells in the given Any. + * + * \param a an Any object. + */ + virtual amrex::IntVect AnyGrowVect (amrex::Any const& a) const override + { + auto const& mfs = a.get(); + return mfs[0].nGrowVect(); + } + + /** + * \brief Copy data from source Any to destination Any. + * + * \param dst destination Any. + * \param src source Any. + * \param ng number of ghost cells included in the operation. + */ + virtual void AnyCopy (amrex::Any& dst, amrex::Any const& src, amrex::IntVect const& ng) const override + { + auto& dmf = dst.get(); + auto const& smf = src.get(); + for (int idim=0; idim < 3; ++idim) { + amrex::MultiFab::Copy(dmf[idim], smf[idim], 0, 0, 1, ng); + } + } + + /** + * \brief Add data from source Any to destination Any. + * + * \param dst destination Any. + * \param src source Any. + * \param ng number of ghost cells included in the operation. + */ + virtual void AnyAdd (amrex::Any& dst, amrex::Any const& src, amrex::IntVect const& ng) const override + { + auto& dmf = dst.get(); + auto const& smf = src.get(); + for (int idim=0; idim < 3; ++idim) { + amrex::MultiFab::Add(dmf[idim], smf[idim], 0, 0, 1, ng); + } + } + + /** + * \brief Set the given Any to zero. + * + * \param a an Any object. + */ + virtual void AnySetToZero (amrex::Any& a) const override + { + auto& mfs = a.get(); + for (int idim=0; idim < 3; ++idim) { + mfs[idim].setVal(amrex::Real(0.0)); + } + } + + /** + * \brief Set boundary (i.e., ghost cells) the given Any to zero. + * + * \param a an Any object. + */ + virtual void AnySetBndryToZero (amrex::Any& a) const override + { + auto& mfs = a.get(); + for (int idim=0; idim < 3; ++idim) { + mfs[idim].setBndry(amrex::Real(0.0), 0, 1); + } + } + +#ifdef AMREX_USE_EB + /** + * \brief Set covered region of the given Any to zero. + * + * \param a an Any object. + */ + virtual void AnySetCoveredToZero (amrex::Any& a) const override + { + auto& mfs = a.get(); + for (int idim=0; idim < 3; ++idim) { + amrex::EB_set_covered(mfs[idim], 0, 1, 0, amrex::Real(0.0)); + } + } +#endif + + /** + * \brief ParallelCopy from source Any ot destination Any. + * + * \param dst destination Any. + * \param src source Any. + * \param src_nghost number of ghost cells in the source included in the operation. + * \param dst_nghost number of ghost cells in the destination included in the operation. + * \param period Periodicity. + */ + virtual void AnyParallelCopy (amrex::Any& dst, amrex::Any const& src, + amrex::IntVect const& src_nghost, amrex::IntVect const& dst_nghost, + amrex::Periodicity const& period = amrex::Periodicity::NonPeriodic()) const override + { + auto& dmf = dst.get(); + auto const& smf = src.get(); + for (int idim=0; idim < 3; ++idim) { + dmf[idim].ParallelCopy_nowait(smf[idim], 0, 0, 1, src_nghost, dst_nghost, period); + } + for (int idim=0; idim < 3; ++idim) { + dmf[idim].ParallelCopy_finish(); + } + } + + /** + * \brief Return the infinity norm of the given Any. + * + * \param a an Any object. + */ + virtual amrex::Real AnyNormInf (amrex::Any& a) const override + { + auto& mfs = a.get(); + amrex::Real r = amrex::Real(0.0); + for (int idim=0; idim < 3; ++idim) { + auto tmp = mfs[idim].norminf(0, 0, true); + r = std::max(r, tmp); + } + amrex::ParallelAllReduce::Max(r, amrex::ParallelContext::CommunicatorSub()); + return r; + } + + /** + * \brief Return the infinity norm of the masked region of the given Any. + * + * For a composite solve with multiple AMR levels, the region covered by + * finer AMR levels are not included in the operation. + * + * \parame amrlev AMR level. + * \param a an Any object. + * \parame local determines if the reduction is local (i.e., no MPI communication) or not. + */ + virtual amrex::Real AnyNormInfMask (int amrlev, amrex::Any const& a, bool local) const override + { + amrex::ignore_unused(amrlev, a, local); + amrex::Abort("TODO: AnyNormInfMask"); + // This is only needed for multi-level composite solve + return amrex::Real(0.0); + } + + /** + * \brief Compute residual of the original form, r = b - Ax. + * + * \param amrlev AMR level + * \param resid residual + * \param x the solution x + * \param b the RHS b + * \param crse_bcdata provides Dirichlet BC at AMR coarse/fine interface. + * It's a nullptr for single level solve. + */ + virtual void AnySolutionResidual (int amrlev, amrex::Any& resid, amrex::Any& x, amrex::Any const& b, + amrex::Any const* crse_bcdata = nullptr) override + { + amrex::ignore_unused(amrlev, resid, x, b, crse_bcdata); + amrex::Abort("TODO: AnySolutionResidual"); + } + + /** + * \brief Compute residual of the residual correction form, r = b - Ax. + * + * \param amrlev AMR level. + * \param resid residual of the residual correction form. + * \param x the correction. + * \param b the RHS for the residual correction form (i.e., the residual of the original form. + * \param bc_mode is either Homogeneous or Inhomogeneous. + * \param crse_bcdata provides inhomogenous Dirichlet BC at AMR coarse/fine interface. + * It's ignored for homogeneous Dirichlet BC. + */ + virtual void AnyCorrectionResidual (int amrlev, int mglev, amrex::Any& resid, amrex::Any& x, + const amrex::Any& b, MLLinOp::BCMode bc_mode, + const amrex::Any* crse_bcdata=nullptr) override + { + amrex::ignore_unused(amrlev, mglev, resid, x, b, bc_mode, crse_bcdata); + amrex::Abort("TODO: AnyCorrectionResidual"); + } + + /** + * \brief Reflux + * + * This modifies the coarse level residual at the coarse/fine interface. + * + * \param crse_amrlev coarse AMR level. + * \param res coarse level residual. + * \param crse_sol coarse level x. + * \param crse_rhs coarse level b. + * \param fine_res fine level residual. This may not be needed depending on the coarse/fine stencil. + * \param fine_sol fine level x. + * \param fine_rhs fine level b. + */ + virtual void AnyReflux (int crse_amrlev, + amrex::Any& res, const amrex::Any& crse_sol, const amrex::Any& crse_rhs, + amrex::Any& fine_res, amrex::Any& fine_sol, const amrex::Any& fine_rhs) override + { + amrex::ignore_unused(crse_amrlev, res, crse_sol, crse_rhs, fine_res, fine_sol, fine_rhs); + amrex::Abort("TODO: AnyReflux"); + // This is only needed for multi-level composite solve + } + + /** + * \brief Average down residual from fine to coarse AMR level. + * + * \param clev coarse ARR level. + * \param cres coarse level residual. + * \param fres fine level residual. + */ + virtual void AnyAvgDownResAmr (int clev, amrex::Any& cres, amrex::Any const& fres) const override + { + amrex::ignore_unused(clev, cres, fres); + amrex::Abort("TODO: AnyAvgDownResAmr"); + // This is only needed for mulit-level composite solve. + // And maybe there is nothing neeed to be done here, like in the nodal projection solver. + } + + /** + * \brief Average down residual from fine to coarse MG level. + * + * This is only needed for MG F-cycle, and we don't need to implement this for V-cycle. + * + * \param clev coarse MG level. + * \param cres coarse level residual. + * \param fres fine level residual. + */ + virtual void AnyAvgDownResMG (int clev, amrex::Any& cres, amrex::Any const& fres) const override + { + amrex::ignore_unused(clev, cres, fres); + amrex::Abort("TODO: AnyAvgDownResMG"); // Not needed for V-cycle. + } + + /** + * \brief Smooth the given level. + * + * \param amrlev AMR level. Note that the lowest level is always 0. + * \param mglev MG level. Note that mglev+1 is one level coarser than mglev. + * \param sol x + * \param rhs b + * \param skip_fillboundary a flag for if we need to fill ghost cells in this function. + */ + virtual void AnySmooth (int amrlev, int mglev, amrex::Any& sol, const amrex::Any& rhs, + bool skip_fillboundary=false) const override + { + amrex::ignore_unused(amrlev, mglev, sol, rhs, skip_fillboundary); + amrex::Abort("TODO: AnySmooth"); + } + + /** + * \brief Restriction from fine to coarse MG level. + * + * \param amrlev AMR level. + * \param cmglev coarse MG level. The fine MG level is cmglev-1. + * \param crse coarse data. + * \param fine fine data. This is not const& because we may need to fill its ghost cells. + */ + virtual void AnyRestriction (int amrlev, int cmglev, amrex::Any& crse, amrex::Any& fine) const override + { + amrex::ignore_unused(amrlev, cmglev, crse, fine); + amrex::Abort("TODO: AnyRestriction"); + } + + /** + * \brief Add interpolated coarse data onto the fine MG level. + * + * Note that it's an ADD operation. + * + * \param amrlev AMR level. + * \param fmglev fine MG level. The coarse MG level is fmglev+1. + * \param fine fine MG level data. + * \param crse coarse MG level data. + */ + virtual void AnyInterpolationMG (int amrlev, int fmglev, amrex::Any& fine, const amrex::Any& crse) const override + { + amrex::ignore_unused(amrlev, fmglev, fine, crse); + amrex::Abort("TODO: AnyInterpolationMG"); + } + + /** + * \brief Assign (i.e., copy) interpolated coarse data onto the fine MG level. + * + * Note that it's an ASSIGN operation. This is used in MG F-cycle, and + * does not need to be implemented for V-cycle. + * + * \param amrlev AMR level. + * \param fmglev fine MG level. The coarse MG level is fmglev+1. + * \param fine fine MG level data. + * \param crse coarse MG level data. + */ + virtual void AnyInterpAssignMG (int amrlev, int fmglev, amrex::Any& fine, amrex::Any& crse) const override + { + amrex::ignore_unused(amrlev, fmglev, fine, crse); + amrex::Abort("TODO: AnyInterpAssignMG"); // not needed for V-cycle. + } + + /** + * \brief Interpolate data from coarse to fine AMR level. + * + * \param famrlev fine AMR level. The coarse AMR level is famrlev-1. + * \param fine data on fine AMR level. + * \param crse data on coarse AMR level. + */ + virtual void AnyInterpolationAmr (int famrlev, amrex::Any& fine, const amrex::Any& crse, + amrex::IntVect const& /*nghost*/) const override + { + amrex::ignore_unused(famrlev, fine, crse); + // This is only needed for multi-level composite solve + amrex::Abort("TODO: AnyInterpolationAmr"); + } + + /** + * \brief Average down x and b from fine to coarse AMR level. + * + * This is called before V-cycle to make data on AMR levels consistent. + * + * \param camrlev coarse AMR level. The fine level is camrlev+1. + * \param crse_sol x on coarse level. + * \param crse_rhs b on coarse level. + * \param fine_sol x on fine level. + * \param fine_rhs b on fine level. + */ + virtual void AnyAverageDownSolutionRHS (int camrlev, amrex::Any& crse_sol, amrex::Any& crse_rhs, + const amrex::Any& fine_sol, const amrex::Any& fine_rhs) override + { + amrex::ignore_unused(camrlev, crse_sol, crse_rhs, fine_sol, fine_rhs); + // This is only needed for multi-level composite solve + amrex::Abort("AnyAverageDownSolutionRHS"); + } + + /** + * \brief Average down and synchronize AMR data. + * + * Synchronize the data on each level. That is the nodal data in the + * same MultiFab needs to be synchronized. This function also needs to + * average down the data from fine to coarse AMR levels. + * + * \param sol data on all AMR levels. + */ + virtual void AnyAverageDownAndSync (amrex::Vector& sol) const override + { + amrex::ignore_unused(sol); + // Even for single level, we shoudl synchronize the data on level 0. + amrex::Abort("TODO: AnyAverageDownAndSync"); + } + + /** + * \brief Prepare the solver for MG cycle. + */ + virtual void prepareForSolve () override + { + amrex::Abort("TODO: prepareForSolve"); + } +}; + +} + + +#endif 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..a1e897e85ba 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; @@ -84,8 +110,6 @@ MLMG::solve (const Vector& a_sol, const Vector& a_rh computeMLResidual(finest_amr_lev); - int ncomp = linop.getNComp(); - bool local = true; Real resnorm0 = MLResNormInf(finest_amr_lev, local); Real rhsnorm0 = MLRhsNormInf(local); @@ -200,9 +224,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], ng_back); } } @@ -229,16 +252,13 @@ 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], nghost); // compute residual for the coarse AMR level computeResWithCrseSolFineCor(alev-1,alev); @@ -250,7 +270,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 +277,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], 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], 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], nghost); } // Update fine AMR level correction @@ -283,14 +305,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], 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], nghost); } } - averageDownAndSync(); + linop.AnyAverageDownAndSync(sol); } // Compute multi-level Residual (res) up to amrlevmax. @@ -301,11 +323,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 +337,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 @@ -333,39 +347,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, 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. @@ -374,20 +377,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, nghost); } void @@ -413,16 +415,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 +433,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 +447,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 +455,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 +465,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 +494,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 +507,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"; } @@ -523,19 +524,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 +537,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], 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], nghost); } } @@ -563,17 +557,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); + IntVect nghost(0); + if (cf_strategy == CFStrategy::ghostnodes) nghost = IntVect(linop.getNGrow(alev)); - const MultiFab& crse_cor = *cor[alev-1][0]; - MultiFab& fine_cor = *cor[alev][0]; - - 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 +572,12 @@ 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, IntVect(ng_dst)); + linop.AnySetToZero(cfine); + linop.AnyParallelCopy(cfine, crse_cor, 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 +588,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) @@ -830,31 +599,24 @@ MLMG::addInterpCorrection (int alev, int mglev) { BL_PROFILE("MLMG::addInterpCorrection()"); - 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, IntVect(0)); + linop.AnyParallelCopy(cfine,crse_cor,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 +627,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 +656,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 +668,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 @@ -914,8 +676,6 @@ MLMG::actualBottomSolve () { BL_PROFILE("MLMG::actualBottomSolve()"); - const int ncomp = linop.getNComp(); - if (!linop.isBottomActive()) return; auto bottom_start_time = amrex::second(); @@ -924,28 +684,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, ng); + linop.AnyCopy(raii_b, b, ng); bottom_b = &raii_b; makeSolvable(amrlev,mglev,*bottom_b); @@ -973,7 +733,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 +743,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 +755,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 +766,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 +787,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 +809,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()"); @@ -1147,7 +828,6 @@ MLMG::prepareForSolve (const Vector& a_sol, 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.AnySetBndryToZero(sol[alev]); + 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, 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], IntVect(0)); + linop.AnySetBndryToZero(sol[alev]); + sol_is_alias[alev] = false; } } @@ -1202,10 +881,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, ng_rhs); } - MultiFab::Copy(rhs[alev], *a_rhs[alev], 0, 0, ncomp, ng_rhs); + linop.AnyCopy(rhs[alev], a_rhs[alev], ng_rhs); linop.applyMetricTerm(alev, 0, rhs[alev]); linop.unimposeNeumannBC(alev, rhs[alev]); linop.applyInhomogNeumannTerm(alev, rhs[alev]); @@ -1215,38 +893,37 @@ 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.AnySetCoveredToZero(rhs[alev]); + linop.AnySetCoveredToZero(sol[alev]); } #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(); } IntVect ng = linop.isCellCentered() ? IntVect(0) : IntVect(1); if (cf_strategy == CFStrategy::ghostnodes) ng = ng_rhs; if (!solve_called) { - linop.make(res, ncomp, ng); - linop.make(rescor, ncomp, ng); + linop.make(res, ng); + linop.make(rescor, ng); } for (int alev = 0; alev <= finest_amr_lev; ++alev) { const int nmglevs = linop.NMGLevels(alev); for (int mglev = 0; mglev < nmglevs; ++mglev) { - res[alev][mglev].setVal(0.0); - rescor[alev][mglev].setVal(0.0); + linop.AnySetToZero(res [alev][mglev]); + linop.AnySetToZero(rescor[alev][mglev]); } } @@ -1261,12 +938,9 @@ MLMG::prepareForSolve (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, _ng); } - cor[alev][mglev]->setVal(0.0); + linop.AnySetToZero(cor[alev][mglev]); } } @@ -1280,12 +954,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, _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 +965,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] = linop.AnyMake(alev, 0, _ng); } - 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)); - } - } - } -#endif + linop.AnySetToZero(cor_hold[alev][0]); } if (linop.m_parent) { @@ -1379,7 +1028,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 +1041,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 +1066,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 +1116,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 +1147,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, 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 +1175,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 +1259,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 +1293,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 = offset.size(); + 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 = offset.size(); for (int c = 0; c < ncomp; ++c) { - mf.plus(-offset[c], c, 1); + amrex::Print() << "MLMG: Subtracting " << offset[c] + << " from mf component c = " << c + << " on level (" << amrlev << ", " << mglev << ")\n"; } -#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"; - } - 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 +1391,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..1935be89f1d 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLNodeLinOp.H @@ -36,10 +36,6 @@ public: const Vector const*>& a_factory = {}, int a_eb_limit_coarsening = -1); - virtual void setLevelBC (int /*amrlev*/, const MultiFab* /*levelbcdata*/, - const MultiFab* = nullptr, const MultiFab* = nullptr, - const MultiFab* = nullptr) final override {} - virtual void apply (int amrlev, int mglev, MultiFab& out, MultiFab& in, BCMode bc_mode, StateMode s_mode, const MLMGBndry* bndry=nullptr) const final override; @@ -59,20 +55,15 @@ 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*/, - const MultiFab* /*crse_bcdata*/=nullptr) final override { - 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 +77,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 +93,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 +143,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