Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add option to use multi-level Hypre #142

Merged
merged 1 commit into from
Aug 18, 2024
Merged

Add option to use multi-level Hypre #142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 15 additions & 2 deletions Projections/hydro_MacProjector.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,10 @@
#include <AMReX_MLEBABecLap.H>
#endif

#ifdef AMREX_USE_HYPRE
#include <AMReX_HypreMLABecLap.H>
#endif

namespace Hydro {

class MacProjector
Expand Down Expand Up @@ -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<amrex::Array<amrex::MultiFab const*,AMREX_SPACEDIM> >& a_beta,
const amrex::Vector<amrex::iMultiFab const*>& a_overset_mask = {});

Expand All @@ -96,7 +100,7 @@ public:
#ifndef AMREX_USE_EB
void initProjector (amrex::Vector<amrex::BoxArray> const& a_grids,
amrex::Vector<amrex::DistributionMapping> 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<amrex::iMultiFab const*>& a_overset_mask = {});

void updateBeta (amrex::Real a_const_beta);
Expand Down Expand Up @@ -183,11 +187,15 @@ private:
amrex::Vector<amrex::Geometry> m_geom;

int m_verbose = 0;
int m_maxiter = 200;

bool m_needs_domain_bcs = true;
amrex::Vector<int> m_needs_level_bcs;
bool m_has_robin = false;

amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> m_lobc;
amrex::Array<amrex::LinOpBCType,AMREX_SPACEDIM> m_hibc;

// Location of umac -- face center vs face centroid
amrex::MLMG::Location m_umac_loc;

Expand All @@ -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<amrex::HypreMLABecLap> m_hypremlabeclap;
#endif
};

}
Expand Down
116 changes: 107 additions & 9 deletions Projections/hydro_MacProjector.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ MacProjector::MacProjector (const Vector<Array<MultiFab*,AMREX_SPACEDIM> >& a_um
}

void MacProjector::initProjector (
const LPInfo& a_lpinfo,
LPInfo a_lpinfo,
const Vector<Array<MultiFab const*,AMREX_SPACEDIM> >& a_beta,
const Vector<iMultiFab const*>& a_overset_mask)
{
Expand All @@ -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) {
Expand Down Expand Up @@ -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<MLABecLaplacian>(m_geom, ba, dm, a_lpinfo);
} else {
Expand All @@ -131,6 +141,13 @@ void MacProjector::initProjector (

m_mlmg = std::make_unique<MLMG>(*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;
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -225,6 +246,11 @@ MacProjector::setDomainBC (const Array<LinOpBCType,AMREX_SPACEDIM>& 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
}


Expand All @@ -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
}


Expand All @@ -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)
{
Expand Down Expand Up @@ -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<BoxArray> grids;
Vector<DistributionMapping> 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<HypreMLABecLap>
(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<Array<MultiFab const*,AMREX_SPACEDIM>> 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] )
{
Expand Down Expand Up @@ -345,8 +426,9 @@ MacProjector::getFluxes (const Vector<Array<MultiFab*,AMREX_SPACEDIM> >& a_flux,
const Vector<MultiFab*>& 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) {
Expand All @@ -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);
Expand All @@ -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 );
Expand All @@ -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);
Expand Down Expand Up @@ -460,7 +545,7 @@ MacProjector::averageDownVelocity ()
#ifndef AMREX_USE_EB
void MacProjector::initProjector (Vector<BoxArray> const& a_grids,
Vector<DistributionMapping> const& a_dmap,
const LPInfo& a_lpinfo, Real const a_const_beta,
LPInfo a_lpinfo, Real const a_const_beta,
const Vector<iMultiFab const*>& a_overset_mask)
{
m_const_beta = a_const_beta;
Expand Down Expand Up @@ -489,6 +574,8 @@ void MacProjector::initProjector (Vector<BoxArray> const& a_grids,
}
}

if (m_use_mlhypre) { a_lpinfo.setMaxCoarseningLevel(0); }

if (a_overset_mask.empty()) {
m_poisson = std::make_unique<MLPoisson>(m_geom, ba, dm, a_lpinfo);
} else {
Expand All @@ -499,6 +586,13 @@ void MacProjector::initProjector (Vector<BoxArray> const& a_grids,

m_mlmg = std::make_unique<MLMG>(*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;
Expand Down Expand Up @@ -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
Expand Down