Multi-level Hypre: Fix bugs (#4089)

The stencil was wrong at the corner where the coarse/fine boundary meets
the domain boundary.

Src/Extern/HYPRE/AMReX_HypreMLABecLap.H

@@ -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,
@@ -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;
@@ -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.

Src/Extern/HYPRE/AMReX_HypreMLABecLap.cpp

@@ -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);
@@ -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;
@@ -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>{};,
@@ -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
@@ -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);
 
@@ -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);
@@ -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)
  {
@@ -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
@@ -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
@@ -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

Src/Extern/HYPRE/AMReX_HypreMLABecLap_2D_K.H

@@ -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
@@ -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);
@@ -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);
@@ -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);
@@ -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);