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

Jump to bottom

Multi-level Hypre: Fix bugs #4089

Merged
merged 3 commits into from
Aug 18, 2024
Merged

Multi-level Hypre: Fix bugs #4089

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e39af9b
Multi-level Hypre: Fix bugs
WeiqunZhang Aug 16, 2024
4badc92
clean up
WeiqunZhang Aug 18, 2024
15eef00
add const
WeiqunZhang Aug 18, 2024
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
3 changes: 3 additions & 0 deletions Src/Extern/HYPRE/AMReX_HypreMLABecLap.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@ public:

void setVerbose (int v) { m_verbose = v; }
void setMaxIter (int v) { m_maxiter = v; }
void setIsSingular (bool v) { m_is_singular = v; }

void setup (Real a_ascalar, Real a_bscalar,
Vector<MultiFab const*> const& a_acoefs,
Expand All @@ -65,6 +66,7 @@ private:

int m_verbose = 0;
int m_maxiter = 200;
bool m_is_singular = false;

Vector<Geometry> m_geom;
Vector<BoxArray> m_grids;
Expand All @@ -87,6 +89,7 @@ private:
Vector<std::unique_ptr<MLMGBndry>> m_bndry;
Vector<std::unique_ptr<BndryRegister>> m_bndry_rhs;
Vector<iMultiFab> m_fine_masks;
Vector<iMultiFab> m_crse_masks;

// For coarse cells at coarse/fine interface. The vector is for AMR
// levels.
Expand Down
57 changes: 36 additions & 21 deletions Src/Extern/HYPRE/AMReX_HypreMLABecLap.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,11 +57,16 @@ HypreMLABecLap::HypreMLABecLap (Vector<Geometry> a_geom,
}

m_fine_masks.resize(m_nlevels-1);
m_crse_masks.resize(m_nlevels-1);
for (int ilev = 0; ilev < m_nlevels-1; ++ilev) {
m_fine_masks[ilev] = amrex::makeFineMask(m_grids[ilev], m_dmap[ilev], IntVect(1),
m_grids[ilev+1], m_ref_ratio[ilev],
m_geom[ilev].periodicity(),
0, 1);
m_crse_masks[ilev].define(m_grids[ilev], m_dmap[ilev], 1, 1);
m_crse_masks[ilev].BuildMask(m_geom[ilev].Domain(),
m_geom[ilev].periodicity(),
1, 0, 0, 1);
}

m_c2f_offset_from.resize(m_nlevels-1);
Expand Down Expand Up @@ -588,12 +593,19 @@ void HypreMLABecLap::setup (Real a_ascalar, Real a_bscalar,
const auto boxlo = amrex::lbound(vbx);
const auto boxhi = amrex::ubound(vbx);
// Set up stencil part of the matrix
auto fixed_pt = IntVect::TheMaxVector();
if (m_is_singular && m_nlevels-1 == ilev) {
auto const& box0 = m_grids.back()[0];
fixed_pt = box0.smallEnd() + 1;
// This cell does not have any non-stencil entries. So it's
// a good point for fixing singularity.
}
amrex::fill(matfab,
[=] AMREX_GPU_HOST_DEVICE (GpuArray<Real,stencil_size>& sten,
int i, int j, int k)
{
hypmlabeclap_mat(sten, i, j, k, boxlo, boxhi, sa, afab, sb, dx, bfabs,
bctype, bcl, bcmsk, bcval, bcrhs, ilev);
bctype, bcl, bcmsk, bcval, bcrhs, ilev, fixed_pt);
});

bool need_sync = true;
Expand Down Expand Up @@ -636,6 +648,7 @@ void HypreMLABecLap::setup (Real a_ascalar, Real a_bscalar,
auto const& c2f_offset_to_a = m_c2f_offset_to[ilev].const_array(mfi);
auto const& mat_a = matfab.array();
auto const& fine_mask = m_fine_masks[ilev].const_array(mfi);
auto const& crse_mask = m_crse_masks[ilev].const_array(mfi);
AMREX_D_TERM(auto offset_bx_a = m_offset_cf_bcoefs[ilev][0].isDefined()
? m_offset_cf_bcoefs[ilev][0].const_array(mfi)
: Array4<int const>{};,
Expand Down Expand Up @@ -664,7 +677,7 @@ void HypreMLABecLap::setup (Real a_ascalar, Real a_bscalar,
c2f_offset_to_a, dx, sb,
AMREX_D_DECL(offset_bx_a,offset_by_a,offset_bz_a),
AMREX_D_DECL(p_bx, p_by, p_bz),
fine_mask,rr);
fine_mask,rr, crse_mask);
});
if (c2f_total_from > 0) {
#ifdef AMREX_USE_GPU
Expand Down Expand Up @@ -838,8 +851,8 @@ void HypreMLABecLap::setup (Real a_ascalar, Real a_bscalar,
HYPRE_SStructSSAMGSetNumPostRelax(m_ss_solver, 4);
HYPRE_SStructSSAMGSetNumCoarseRelax(m_ss_solver, 4);

HYPRE_SStructSSAMGSetLogging(m_ss_solver, m_verbose);
HYPRE_SStructSSAMGSetPrintLevel(m_ss_solver, m_verbose);
HYPRE_SStructSSAMGSetLogging(m_ss_solver, 1);
// HYPRE_SStructSSAMGSetPrintLevel(m_ss_solver, 1); /* 0: no, 1: setup, 2: solve, 3:both

// HYPRE_SStructSSAMGSetup(m_ss_solver, A, b, x);

Expand All @@ -854,15 +867,15 @@ void HypreMLABecLap::setup (Real a_ascalar, Real a_bscalar,
HYPRE_BoomerAMGCreate(&m_solver);

HYPRE_BoomerAMGSetOldDefault(m_solver); // Falgout coarsening with modified classical interpolation
HYPRE_BoomerAMGSetStrongThreshold(m_solver, (AMREX_SPACEDIM == 3) ? 0.6 : 0.25); // default is 0.25
HYPRE_BoomerAMGSetStrongThreshold(m_solver, (AMREX_SPACEDIM == 3) ? 0.4 : 0.25); // default is 0.25
HYPRE_BoomerAMGSetRelaxOrder(m_solver, 1); /* 0: default, natural order, 1: C/F relaxation order */
HYPRE_BoomerAMGSetNumSweeps(m_solver, 2); /* Sweeps on fine levels */
// HYPRE_BoomerAMGSetFCycle(m_solver, 1); // default is 0
// HYPRE_BoomerAMGSetCoarsenType(m_solver, 6);
// HYPRE_BoomerAMGSetRelaxType(m_solver, 6); /* G-S/Jacobi hybrid relaxation */

HYPRE_BoomerAMGSetLogging(m_solver, m_verbose);
HYPRE_BoomerAMGSetPrintLevel(m_solver, m_verbose);
HYPRE_BoomerAMGSetLogging(m_solver, 1);
// HYPRE_BoomerAMGSetPrintLevel(m_solver, 1); /* 0: no, 1: setup, 2: solve, 3:both

HYPRE_ParCSRMatrix par_A;
HYPRE_SStructMatrixGetObject(m_ss_A, (void**) &par_A);
Expand Down Expand Up @@ -956,6 +969,9 @@ void HypreMLABecLap::solve (Vector<MultiFab*> const& a_sol, Vector<MultiFab cons
amrex::Abort("HypreMLABecLap::solve: TODO abstol > 0");
}

HYPRE_Int num_iterations;
Real final_res;

#ifdef AMREX_FEATURE_HYPRE_SSAMG
if (m_hypre_solver_id == HypreSolverID::SSAMG)
{
Expand All @@ -965,15 +981,13 @@ void HypreMLABecLap::solve (Vector<MultiFab*> const& a_sol, Vector<MultiFab cons

HYPRE_SStructSSAMGSolve(m_ss_solver, m_ss_A, m_ss_b, m_ss_x);

if (m_verbose) {
HYPRE_Int num_iterations;
Real res;
HYPRE_SStructSSAMGGetNumIterations(m_ss_solver, &num_iterations);
HYPRE_SStructSSAMGGetFinalRelativeResidualNorm(m_ss_solver, &res);
HYPRE_SStructSSAMGGetNumIterations(m_ss_solver, &num_iterations);
HYPRE_SStructSSAMGGetFinalRelativeResidualNorm(m_ss_solver, &final_res);

if (m_verbose) {
amrex::Print() << "\n" << num_iterations
<< " Hypre SSAMG Iterations, Relative Residual "
<< res << '\n';
<< final_res << '\n';
}
} else
#endif
Expand All @@ -990,17 +1004,20 @@ void HypreMLABecLap::solve (Vector<MultiFab*> const& a_sol, Vector<MultiFab cons

HYPRE_BoomerAMGSolve(m_solver, par_A, par_b, par_x);

if (m_verbose) {
HYPRE_Int num_iterations;
Real res;
HYPRE_BoomerAMGGetNumIterations(m_solver, &num_iterations);
HYPRE_BoomerAMGGetFinalRelativeResidualNorm(m_solver, &res);
HYPRE_BoomerAMGGetNumIterations(m_solver, &num_iterations);
HYPRE_BoomerAMGGetFinalRelativeResidualNorm(m_solver, &final_res);

if (m_verbose) {
amrex::Print() << "\n" << num_iterations
<< " Hypre SS BoomerAMG Iterations, Relative Residual "
<< res << '\n';
<< final_res << '\n';
}
}

if (final_res > reltol) {
amrex::Abort("Hypre failed to converge after "+std::to_string(num_iterations)+
" iterations. Final relative residual is "+std::to_string(final_res));
}
}

// Get solution
Expand Down Expand Up @@ -1044,8 +1061,6 @@ void HypreMLABecLap::solve (Vector<MultiFab*> const& a_sol, Vector<MultiFab cons
for (int ilev = m_nlevels-2; ilev >= 0; --ilev) {
amrex::average_down(*a_sol[ilev+1], *a_sol[ilev], 0, ncomp, m_ref_ratio[ilev]);
}

// xxxxx abort if convergence is not reached.
}

#ifdef AMREX_USE_GPU
Expand Down
34 changes: 25 additions & 9 deletions Src/Extern/HYPRE/AMReX_HypreMLABecLap_2D_K.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ void hypmlabeclap_c2f (int i, int j, int k,
Array4<int const> const& offset_by,
Real const* bx, Real const* by,
Array4<int const> const& fine_mask,
IntVect const& rr)
IntVect const& rr, Array4<int const> const& crse_mask)
{
if (fine_mask(i,j,k)) {
// Let's set off-diagonal elements to zero
Expand Down Expand Up @@ -159,9 +159,13 @@ void hypmlabeclap_c2f (int i, int j, int k,
Real xInt = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
Real xc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
if (fine_mask(i-1,j,k)) {
if ( fine_mask(i-1,j,k) ||
!crse_mask(i-1,j,k))
{
poly_interp_coeff<2>(xInt, &(xc[1]), &(ct[1]));
} else if (fine_mask(i+1,j,k)) {
} else if ( fine_mask(i+1,j,k) ||
!crse_mask(i+1,j,k))
{
poly_interp_coeff<2>(xInt, xc, ct);
} else {
poly_interp_coeff<3>(xInt, xc, ct);
Expand Down Expand Up @@ -202,9 +206,13 @@ void hypmlabeclap_c2f (int i, int j, int k,
Real yInt = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
Real yc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
if (fine_mask(i,j-1,k)) {
if ( fine_mask(i,j-1,k) ||
!crse_mask(i,j-1,k))
{
poly_interp_coeff<2>(yInt, &(yc[1]), &(ct[1]));
} else if (fine_mask(i,j+1,k)) {
} else if ( fine_mask(i,j+1,k) ||
!crse_mask(i,j+1,k))
{
poly_interp_coeff<2>(yInt, yc, ct);
} else {
poly_interp_coeff<3>(yInt, yc, ct);
Expand Down Expand Up @@ -244,9 +252,13 @@ void hypmlabeclap_c2f (int i, int j, int k,
Real yInt = Real(-0.5) + (iry+Real(0.5))/Real(rr[1]);
Real yc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
if (fine_mask(i,j-1,k)) {
if ( fine_mask(i,j-1,k) ||
!crse_mask(i,j-1,k))
{
poly_interp_coeff<2>(yInt, &(yc[1]), &(ct[1]));
} else if (fine_mask(i,j+1,k)) {
} else if ( fine_mask(i,j+1,k) ||
!crse_mask(i,j+1,k))
{
poly_interp_coeff<2>(yInt, yc, ct);
} else {
poly_interp_coeff<3>(yInt, yc, ct);
Expand Down Expand Up @@ -286,9 +298,13 @@ void hypmlabeclap_c2f (int i, int j, int k,
Real xInt = Real(-0.5) + (irx+Real(0.5))/Real(rr[0]);
Real xc[3] = {Real(-1.0), Real(0.0), Real(1.0)};
Real ct[3] = {Real(0.0), Real(0.0), Real(0.0)};
if (fine_mask(i-1,j,k)) {
if ( fine_mask(i-1,j,k) ||
!crse_mask(i-1,j,k))
{
poly_interp_coeff<2>(xInt, &(xc[1]), &(ct[1]));
} else if (fine_mask(i+1,j,k)) {
} else if ( fine_mask(i+1,j,k) ||
!crse_mask(i+1,j,k))
{
poly_interp_coeff<2>(xInt, xc, ct);
} else {
poly_interp_coeff<3>(xInt, xc, ct);
Expand Down
Loading
Loading