diff --git a/Projections/hydro_MacProjector.H b/Projections/hydro_MacProjector.H index df6c7bf85..85301c6e5 100644 --- a/Projections/hydro_MacProjector.H +++ b/Projections/hydro_MacProjector.H @@ -10,6 +10,10 @@ #include #endif +#ifdef AMREX_USE_HYPRE +#include +#endif + namespace Hydro { class MacProjector @@ -84,7 +88,7 @@ public: /** Initialize the underlying linear operator and MLMG instances */ void initProjector ( - const amrex::LPInfo& a_lpinfo, + amrex::LPInfo a_lpinfo, const amrex::Vector >& a_beta, const amrex::Vector& a_overset_mask = {}); @@ -96,7 +100,7 @@ public: #ifndef AMREX_USE_EB void initProjector (amrex::Vector const& a_grids, amrex::Vector const& a_dmap, - const amrex::LPInfo& a_lpinfo, amrex::Real a_const_beta, + amrex::LPInfo a_lpinfo, amrex::Real a_const_beta, const amrex::Vector& a_overset_mask = {}); void updateBeta (amrex::Real a_const_beta); @@ -183,11 +187,15 @@ private: amrex::Vector m_geom; int m_verbose = 0; + int m_maxiter = 200; bool m_needs_domain_bcs = true; amrex::Vector m_needs_level_bcs; bool m_has_robin = false; + amrex::Array m_lobc; + amrex::Array m_hibc; + // Location of umac -- face center vs face centroid amrex::MLMG::Location m_umac_loc; @@ -199,6 +207,11 @@ private: amrex::MLMG::Location m_divu_loc; bool m_needs_init = true; + + bool m_use_mlhypre = false; +#ifdef AMREX_USE_HYPRE + std::unique_ptr m_hypremlabeclap; +#endif }; } diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp index 8dd8fafc3..c60e0a130 100644 --- a/Projections/hydro_MacProjector.cpp +++ b/Projections/hydro_MacProjector.cpp @@ -51,7 +51,7 @@ MacProjector::MacProjector (const Vector >& a_um } void MacProjector::initProjector ( - const LPInfo& a_lpinfo, + LPInfo a_lpinfo, const Vector >& a_beta, const Vector& a_overset_mask) { @@ -69,9 +69,17 @@ void MacProjector::initProjector ( m_fluxes.resize(nlevs); m_divu.resize(nlevs); +#ifdef AMREX_USE_HYPRE + { + ParmParse pp("mac_proj"); + pp.query("use_mlhypre", m_use_mlhypre); + } +#endif + #ifdef AMREX_USE_EB bool has_eb = a_beta[0][0]->hasEBFabFactory(); if (has_eb) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false == m_use_mlhypre, "mlhypre does not work with EB"); m_eb_vel.resize(nlevs); m_eb_factory.resize(nlevs, nullptr); for (int ilev = 0; ilev < nlevs; ++ilev) { @@ -115,6 +123,8 @@ void MacProjector::initProjector ( } } + if (m_use_mlhypre) { a_lpinfo.setMaxCoarseningLevel(0); } + if (a_overset_mask.empty()) { m_abeclap = std::make_unique(m_geom, ba, dm, a_lpinfo); } else { @@ -131,6 +141,13 @@ void MacProjector::initProjector ( m_mlmg = std::make_unique(*m_linop); +#ifdef AMREX_USE_HYPRE + if (m_use_mlhypre) { + m_abeclap->setInterpBndryHalfWidth(1); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_overset_mask.empty(), "mlhypre does not support overset mask yet"); + } +#endif + setOptions(); m_needs_init = false; @@ -183,6 +200,10 @@ void MacProjector::updateCoeffs ( for (int ilev=0; ilev < nlevs; ++ilev) m_abeclap->setBCoeffs(ilev, a_beta[ilev]); } + +#ifdef AMREX_USE_HYPRE + m_hypremlabeclap.reset(); +#endif } void MacProjector::setUMAC( @@ -225,6 +246,11 @@ MacProjector::setDomainBC (const Array& lobc, "MacProjector::setDomainBC: initProjector must be called before calling this method"); m_linop->setDomainBC(lobc, hibc); m_needs_domain_bcs = false; + m_lobc = lobc; + m_hibc = hibc; +#ifdef AMREX_USE_HYPRE + m_hypremlabeclap.reset(); +#endif } @@ -237,6 +263,9 @@ MacProjector::setLevelBC (int amrlev, const MultiFab* levelbcdata, const MultiFa m_linop->setLevelBC(amrlev, levelbcdata, robin_a, robin_b, robin_f); m_needs_level_bcs[amrlev] = false; if (robin_a) { m_has_robin = true; } +#ifdef AMREX_USE_HYPRE + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false == m_use_mlhypre, "mlhypre does not support setLevelBC"); +#endif } @@ -252,8 +281,9 @@ MacProjector::project (Real reltol, Real atol) } } - if ( m_umac[0][0] ) - averageDownVelocity(); + if ( m_umac[0][0] ) { + averageDownVelocity(); + } for (int ilev = 0; ilev < nlevs; ++ilev) { @@ -301,7 +331,58 @@ MacProjector::project (Real reltol, Real atol) m_phi[ilev].setVal(0.0); } - m_mlmg->solve(amrex::GetVecOfPtrs(m_phi), amrex::GetVecOfConstPtrs(m_rhs), reltol, atol); +#ifdef AMREX_USE_HYPRE + if (m_use_mlhypre) { + // We use mlmg to compute the initial residual. It also makes it + // ready for getting fluxes. + m_mlmg->solve(amrex::GetVecOfPtrs(m_phi), amrex::GetVecOfConstPtrs(m_rhs), 1.e10, 0.0); + auto resnorm0 = m_mlmg->getInitResidual(); + + if (m_verbose) { + amrex::Print() << "Initial residual: " << resnorm0 << "\n"; + } + + if (resnorm0 > atol && resnorm0 > Real(0.0)) { + if (!m_hypremlabeclap) { + Vector grids; + Vector dmap; + grids.reserve(m_geom.size()); + dmap.reserve(m_geom.size()); + for (auto const& mf : m_rhs) { + grids.push_back(mf.boxArray()); + dmap.push_back(mf.DistributionMap()); + } + m_hypremlabeclap = std::make_unique + (m_geom, grids, dmap, HypreSolverID::BoomerAMG); + m_hypremlabeclap->setVerbose(m_verbose); + m_hypremlabeclap->setMaxIter(m_maxiter); + m_hypremlabeclap->setIsSingular(m_linop->isSingular(0)); + if (m_poisson) { + m_hypremlabeclap->setup(Real(0.0), Real(-1.0), {}, {}, m_lobc, m_hibc, + amrex::GetVecOfConstPtrs(m_phi)); + } else { + Vector> bcoefs(m_geom.size()); + for (int ilev = 0; ilev < nlevs; ++ilev) { + bcoefs[ilev] = m_abeclap->getBCoeffs(ilev,0); + } + m_hypremlabeclap->setup(Real(0.0), Real(1.0), {}, bcoefs, m_lobc, m_hibc, + amrex::GetVecOfConstPtrs(m_phi)); + } + } + + reltol = std::max(reltol, atol / resnorm0); + atol = 0; + m_hypremlabeclap->solve(amrex::GetVecOfPtrs(m_phi), + amrex::GetVecOfConstPtrs(m_rhs), + reltol, atol); + // Need to prepare mlmg for getFluxes + m_mlmg->prepareForFluxes(amrex::GetVecOfConstPtrs(m_phi)); + } + } else +#endif + { + m_mlmg->solve(amrex::GetVecOfPtrs(m_phi), amrex::GetVecOfConstPtrs(m_rhs), reltol, atol); + } if ( m_umac[0][0] ) { @@ -345,8 +426,9 @@ MacProjector::getFluxes (const Vector >& a_flux, const Vector& a_sol, MLMG::Location a_loc) const { int ilev = 0; - if (m_needs_level_bcs[ilev]) + if (m_needs_level_bcs[ilev]) { m_linop->setLevelBC(ilev, nullptr); + } m_linop->getFluxes(a_flux, a_sol, a_loc); if (m_poisson) { @@ -367,7 +449,6 @@ MacProjector::setOptions () // Default values int maxorder(3); int bottom_verbose(0); - int maxiter(200); int bottom_maxiter(200); Real bottom_rtol(1.0e-4_rt); Real bottom_atol(-1.0_rt); @@ -382,7 +463,7 @@ MacProjector::setOptions () pp.query( "verbose" , m_verbose ); pp.query( "maxorder" , maxorder ); pp.query( "bottom_verbose", bottom_verbose ); - pp.query( "maxiter" , maxiter ); + pp.query( "maxiter" , m_maxiter ); pp.query( "bottom_maxiter", bottom_maxiter ); pp.query( "bottom_rtol" , bottom_rtol ); pp.query( "bottom_atol" , bottom_atol ); @@ -392,11 +473,15 @@ MacProjector::setOptions () pp.query( "num_post_smooth" , num_post_smooth ); pp.query( "num_final_smooth" , num_final_smooth ); + if (m_use_mlhypre) { + maxorder = 3; + } + // Set default/input values m_linop->setMaxOrder(maxorder); m_mlmg->setVerbose(m_verbose); m_mlmg->setBottomVerbose(bottom_verbose); - m_mlmg->setMaxIter(maxiter); + m_mlmg->setMaxIter(m_maxiter); m_mlmg->setBottomMaxIter(bottom_maxiter); m_mlmg->setBottomTolerance(bottom_rtol); m_mlmg->setBottomToleranceAbs(bottom_atol); @@ -460,7 +545,7 @@ MacProjector::averageDownVelocity () #ifndef AMREX_USE_EB void MacProjector::initProjector (Vector const& a_grids, Vector const& a_dmap, - const LPInfo& a_lpinfo, Real const a_const_beta, + LPInfo a_lpinfo, Real const a_const_beta, const Vector& a_overset_mask) { m_const_beta = a_const_beta; @@ -489,6 +574,8 @@ void MacProjector::initProjector (Vector const& a_grids, } } + if (m_use_mlhypre) { a_lpinfo.setMaxCoarseningLevel(0); } + if (a_overset_mask.empty()) { m_poisson = std::make_unique(m_geom, ba, dm, a_lpinfo); } else { @@ -499,6 +586,13 @@ void MacProjector::initProjector (Vector const& a_grids, m_mlmg = std::make_unique(*m_linop); +#ifdef AMREX_USE_HYPRE + if (m_use_mlhypre) { + m_poisson->setInterpBndryHalfWidth(1); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(a_overset_mask.empty(), "mlhypre does not support overset mask yet"); + } +#endif + setOptions(); m_needs_init = false; @@ -540,6 +634,10 @@ void MacProjector::updateBeta (Real a_const_beta) "MacProjector::updateBeta: should not be called for variable beta"); m_const_beta = a_const_beta; + +#ifdef AMREX_USE_HYPRE + m_hypremlabeclap.reset(); +#endif } #endif