From 8222a13de06d3fe063a1f635623348905c2475f8 Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Sat, 17 Aug 2024 21:11:00 -0500 Subject: [PATCH] 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 | 3 + Src/Extern/HYPRE/AMReX_HypreMLABecLap.cpp | 57 +++-- Src/Extern/HYPRE/AMReX_HypreMLABecLap_2D_K.H | 34 ++- Src/Extern/HYPRE/AMReX_HypreMLABecLap_3D_K.H | 224 +++++++++++-------- Src/Extern/HYPRE/AMReX_HypreMLABecLap_K.H | 8 +- Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H | 33 ++- Src/LinearSolvers/MLMG/AMReX_MLLinOp.H | 2 + Src/LinearSolvers/MLMG/AMReX_MLMG.H | 12 + 8 files changed, 250 insertions(+), 123 deletions(-) diff --git a/Src/Extern/HYPRE/AMReX_HypreMLABecLap.H b/Src/Extern/HYPRE/AMReX_HypreMLABecLap.H index 04147207d36..6f687766b86 100644 --- a/Src/Extern/HYPRE/AMReX_HypreMLABecLap.H +++ b/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 const& a_acoefs, @@ -65,6 +66,7 @@ private: int m_verbose = 0; int m_maxiter = 200; + bool m_is_singular = false; Vector m_geom; Vector m_grids; @@ -87,6 +89,7 @@ private: Vector> m_bndry; Vector> m_bndry_rhs; Vector m_fine_masks; + Vector m_crse_masks; // For coarse cells at coarse/fine interface. The vector is for AMR // levels. diff --git a/Src/Extern/HYPRE/AMReX_HypreMLABecLap.cpp b/Src/Extern/HYPRE/AMReX_HypreMLABecLap.cpp index d7621c6bd0e..b726f307bba 100644 --- a/Src/Extern/HYPRE/AMReX_HypreMLABecLap.cpp +++ b/Src/Extern/HYPRE/AMReX_HypreMLABecLap.cpp @@ -57,11 +57,16 @@ HypreMLABecLap::HypreMLABecLap (Vector 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& 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{};, @@ -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 const& a_sol, Vector 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 const& a_sol, Vector const& a_sol, Vector 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 const& a_sol, Vector= 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 diff --git a/Src/Extern/HYPRE/AMReX_HypreMLABecLap_2D_K.H b/Src/Extern/HYPRE/AMReX_HypreMLABecLap_2D_K.H index 57a37f19bf3..7d083e7d98c 100644 --- a/Src/Extern/HYPRE/AMReX_HypreMLABecLap_2D_K.H +++ b/Src/Extern/HYPRE/AMReX_HypreMLABecLap_2D_K.H @@ -109,7 +109,7 @@ void hypmlabeclap_c2f (int i, int j, int k, Array4 const& offset_by, Real const* bx, Real const* by, Array4 const& fine_mask, - IntVect const& rr) + IntVect const& rr, Array4 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); diff --git a/Src/Extern/HYPRE/AMReX_HypreMLABecLap_3D_K.H b/Src/Extern/HYPRE/AMReX_HypreMLABecLap_3D_K.H index 8e6e1a39b14..431650236f7 100644 --- a/Src/Extern/HYPRE/AMReX_HypreMLABecLap_3D_K.H +++ b/Src/Extern/HYPRE/AMReX_HypreMLABecLap_3D_K.H @@ -166,7 +166,7 @@ void hypmlabeclap_c2f (int i, int j, int k, Array4 const& offset_bz, Real const* bx, Real const* by, Real const* bz, Array4 const& fine_mask, - IntVect const& rr) + IntVect const& rr, Array4 const& crse_mask) { if (fine_mask(i,j,k)) { // Let's set off-diagonal elements to zero @@ -191,7 +191,11 @@ void hypmlabeclap_c2f (int i, int j, int k, (! fine_mask(i,j-1,k-1)) && (! fine_mask(i,j+1,k-1)) && (! fine_mask(i,j-1,k+1)) && - (! fine_mask(i,j+1,k+1))) + (! fine_mask(i,j+1,k+1)) && + ( crse_mask(i,j-1,k-1)) && + ( crse_mask(i,j+1,k-1)) && + ( crse_mask(i,j-1,k+1)) && + ( crse_mask(i,j+1,k+1))) { corner[0] = true; } @@ -199,7 +203,11 @@ void hypmlabeclap_c2f (int i, int j, int k, (! fine_mask(i-1,j,k-1)) && (! fine_mask(i+1,j,k-1)) && (! fine_mask(i-1,j,k+1)) && - (! fine_mask(i+1,j,k+1))) + (! fine_mask(i+1,j,k+1)) && + ( crse_mask(i-1,j,k-1)) && + ( crse_mask(i+1,j,k-1)) && + ( crse_mask(i-1,j,k+1)) && + ( crse_mask(i+1,j,k+1))) { corner[1] = true; } @@ -207,7 +215,11 @@ void hypmlabeclap_c2f (int i, int j, int k, (! fine_mask(i-1,j-1,k)) && (! fine_mask(i+1,j-1,k)) && (! fine_mask(i-1,j+1,k)) && - (! fine_mask(i+1,j+1,k))) + (! fine_mask(i+1,j+1,k)) && + ( crse_mask(i-1,j-1,k)) && + ( crse_mask(i+1,j-1,k)) && + ( crse_mask(i-1,j+1,k)) && + ( crse_mask(i+1,j+1,k))) { corner[2] = true; } @@ -253,28 +265,34 @@ void hypmlabeclap_c2f (int i, int j, int k, Real fac0 = fac*cc[0]; Real s0 = Real(1.0); - if (!fine_mask(i-1,j,k) && !fine_mask(i+1,j,k)) { - s0 -= x*x; - stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); - stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); - } else if (!fine_mask(i-1,j,k)) { + if ( fine_mask(i-1,j,k) || + !crse_mask(i-1,j,k)) + { + s0 += Real(-0.5)*x; + stencil(i,j,k)[2] += fac0*Real(0.5)*x; + } else if ( fine_mask(i+1,j,k) || + !crse_mask(i+1,j,k)) { s0 += Real(0.5)*x; stencil(i,j,k)[1] += fac0*Real(-0.5)*x; } else { - s0 += Real(-0.5)*x; - stencil(i,j,k)[2] += fac0*Real(0.5)*x; + s0 -= x*x; + stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); + stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); } - if (!fine_mask(i,j-1,k) && !fine_mask(i,j+1,k)) { - s0 -= y*y; - stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); - stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); - } else if (!fine_mask(i,j-1,k)) { + if ( fine_mask(i,j-1,k) || + !crse_mask(i,j-1,k)) + { + s0 += Real(-0.5)*y; + stencil(i,j,k)[4] += fac0*Real(0.5)*y; + } else if ( fine_mask(i,j+1,k) || + !crse_mask(i,j+1,k)) { s0 += Real(0.5)*y; stencil(i,j,k)[3] += fac0*Real(-0.5)*y; } else { - s0 += Real(-0.5)*y; - stencil(i,j,k)[4] += fac0*Real(0.5)*y; + s0 -= y*y; + stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); + stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); } stencil(i,j,k)[0] += fac0*s0; @@ -322,28 +340,34 @@ void hypmlabeclap_c2f (int i, int j, int k, Real fac0 = fac*cc[0]; Real s0 = Real(1.0); - if (!fine_mask(i-1,j,k) && !fine_mask(i+1,j,k)) { - s0 -= x*x; - stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); - stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); - } else if (!fine_mask(i-1,j,k)) { + if ( fine_mask(i-1,j,k) || + !crse_mask(i-1,j,k)) + { + s0 += Real(-0.5)*x; + stencil(i,j,k)[2] += fac0*Real(0.5)*x; + } else if ( fine_mask(i+1,j,k) || + !crse_mask(i+1,j,k)) { s0 += Real(0.5)*x; stencil(i,j,k)[1] += fac0*Real(-0.5)*x; } else { - s0 += Real(-0.5)*x; - stencil(i,j,k)[2] += fac0*Real(0.5)*x; + s0 -= x*x; + stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); + stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); } - if (!fine_mask(i,j,k-1) && !fine_mask(i,j,k+1)) { - s0 -= z*z; - stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); - stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); - } else if (!fine_mask(i,j,k-1)) { + if ( fine_mask(i,j,k-1) || + !crse_mask(i,j,k-1)) + { + s0 += Real(-0.5)*z; + stencil(i,j,k)[6] += fac0*Real(0.5)*z; + } else if ( fine_mask(i,j,k+1) || + !crse_mask(i,j,k+1)) { s0 += Real(0.5)*z; stencil(i,j,k)[5] += fac0*Real(-0.5)*z; } else { - s0 += Real(-0.5)*z; - stencil(i,j,k)[6] += fac0*Real(0.5)*z; + s0 -= z*z; + stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); + stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); } stencil(i,j,k)[0] += fac0*s0; @@ -393,28 +417,34 @@ void hypmlabeclap_c2f (int i, int j, int k, Real fac0 = fac*cc[0]; Real s0 = Real(1.0); - if (!fine_mask(i,j-1,k) && !fine_mask(i,j+1,k)) { - s0 -= y*y; - stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); - stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); - } else if (!fine_mask(i,j-1,k)) { + if ( fine_mask(i,j-1,k) || + !crse_mask(i,j-1,k)) + { + s0 += Real(-0.5)*y; + stencil(i,j,k)[4] += fac0*Real(0.5)*y; + } else if ( fine_mask(i,j+1,k) || + !crse_mask(i,j+1,k)) { s0 += Real(0.5)*y; stencil(i,j,k)[3] += fac0*Real(-0.5)*y; } else { - s0 += Real(-0.5)*y; - stencil(i,j,k)[4] += fac0*Real(0.5)*y; + s0 -= y*y; + stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); + stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); } - if (!fine_mask(i,j,k-1) && !fine_mask(i,j,k+1)) { - s0 -= z*z; - stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); - stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); - } else if (!fine_mask(i,j,k-1)) { + if ( fine_mask(i,j,k-1) || + !crse_mask(i,j,k-1)) + { + s0 += Real(-0.5)*z; + stencil(i,j,k)[6] += fac0*Real(0.5)*z; + } else if ( fine_mask(i,j,k+1) || + !crse_mask(i,j,k+1)) { s0 += Real(0.5)*z; stencil(i,j,k)[5] += fac0*Real(-0.5)*z; } else { - s0 += Real(-0.5)*z; - stencil(i,j,k)[6] += fac0*Real(0.5)*z; + s0 -= z*z; + stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); + stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); } stencil(i,j,k)[0] += fac0*s0; @@ -463,28 +493,34 @@ void hypmlabeclap_c2f (int i, int j, int k, Real fac0 = fac*cc[0]; Real s0 = Real(1.0); - if (!fine_mask(i,j-1,k) && !fine_mask(i,j+1,k)) { - s0 -= y*y; - stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); - stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); - } else if (!fine_mask(i,j-1,k)) { + if ( fine_mask(i,j-1,k) || + !crse_mask(i,j-1,k)) + { + s0 += Real(-0.5)*y; + stencil(i,j,k)[4] += fac0*Real(0.5)*y; + } else if ( fine_mask(i,j+1,k) || + !crse_mask(i,j+1,k)) { s0 += Real(0.5)*y; stencil(i,j,k)[3] += fac0*Real(-0.5)*y; } else { - s0 += Real(-0.5)*y; - stencil(i,j,k)[4] += fac0*Real(0.5)*y; + s0 -= y*y; + stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); + stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); } - if (!fine_mask(i,j,k-1) && !fine_mask(i,j,k+1)) { - s0 -= z*z; - stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); - stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); - } else if (!fine_mask(i,j,k-1)) { + if ( fine_mask(i,j,k-1) || + !crse_mask(i,j,k-1)) + { + s0 += Real(-0.5)*z; + stencil(i,j,k)[6] += fac0*Real(0.5)*z; + } else if ( fine_mask(i,j,k+1) || + !crse_mask(i,j,k+1)) { s0 += Real(0.5)*z; stencil(i,j,k)[5] += fac0*Real(-0.5)*z; } else { - s0 += Real(-0.5)*z; - stencil(i,j,k)[6] += fac0*Real(0.5)*z; + s0 -= z*z; + stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); + stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); } stencil(i,j,k)[0] += fac0*s0; @@ -534,28 +570,34 @@ void hypmlabeclap_c2f (int i, int j, int k, Real fac0 = fac*cc[0]; Real s0 = Real(1.0); - if (!fine_mask(i-1,j,k) && !fine_mask(i+1,j,k)) { - s0 -= x*x; - stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); - stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); - } else if (!fine_mask(i-1,j,k)) { + if ( fine_mask(i-1,j,k) || + !crse_mask(i-1,j,k)) + { + s0 += Real(-0.5)*x; + stencil(i,j,k)[2] += fac0*Real(0.5)*x; + } else if ( fine_mask(i+1,j,k) || + !crse_mask(i+1,j,k)) { s0 += Real(0.5)*x; stencil(i,j,k)[1] += fac0*Real(-0.5)*x; } else { - s0 += Real(-0.5)*x; - stencil(i,j,k)[2] += fac0*Real(0.5)*x; + s0 -= x*x; + stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); + stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); } - if (!fine_mask(i,j,k-1) && !fine_mask(i,j,k+1)) { - s0 -= z*z; - stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); - stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); - } else if (!fine_mask(i,j,k-1)) { + if ( fine_mask(i,j,k-1) || + !crse_mask(i,j,k-1)) + { + s0 += Real(-0.5)*z; + stencil(i,j,k)[6] += fac0*Real(0.5)*z; + } else if ( fine_mask(i,j,k+1) || + !crse_mask(i,j,k+1)) { s0 += Real(0.5)*z; stencil(i,j,k)[5] += fac0*Real(-0.5)*z; } else { - s0 += Real(-0.5)*z; - stencil(i,j,k)[6] += fac0*Real(0.5)*z; + s0 -= z*z; + stencil(i,j,k)[5] += fac0*Real(0.5)*z*(z-Real(1.0)); + stencil(i,j,k)[6] += fac0*Real(0.5)*z*(z+Real(1.0)); } stencil(i,j,k)[0] += fac0*s0; @@ -605,28 +647,34 @@ void hypmlabeclap_c2f (int i, int j, int k, Real fac0 = fac*cc[0]; Real s0 = Real(1.0); - if (!fine_mask(i-1,j,k) && !fine_mask(i+1,j,k)) { - s0 -= x*x; - stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); - stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); - } else if (!fine_mask(i-1,j,k)) { + if ( fine_mask(i-1,j,k) || + !crse_mask(i-1,j,k)) + { + s0 += Real(-0.5)*x; + stencil(i,j,k)[2] += fac0*Real(0.5)*x; + } else if ( fine_mask(i+1,j,k) || + !crse_mask(i+1,j,k)) { s0 += Real(0.5)*x; stencil(i,j,k)[1] += fac0*Real(-0.5)*x; } else { - s0 += Real(-0.5)*x; - stencil(i,j,k)[2] += fac0*Real(0.5)*x; + s0 -= x*x; + stencil(i,j,k)[1] += fac0*Real(0.5)*x*(x-Real(1.0)); + stencil(i,j,k)[2] += fac0*Real(0.5)*x*(x+Real(1.0)); } - if (!fine_mask(i,j-1,k) && !fine_mask(i,j+1,k)) { - s0 -= y*y; - stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); - stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); - } else if (!fine_mask(i,j-1,k)) { + if ( fine_mask(i,j-1,k) || + !crse_mask(i,j-1,k)) + { + s0 += Real(-0.5)*y; + stencil(i,j,k)[4] += fac0*Real(0.5)*y; + } else if ( fine_mask(i,j+1,k) || + !crse_mask(i,j+1,k)) { s0 += Real(0.5)*y; stencil(i,j,k)[3] += fac0*Real(-0.5)*y; } else { - s0 += Real(-0.5)*y; - stencil(i,j,k)[4] += fac0*Real(0.5)*y; + s0 -= y*y; + stencil(i,j,k)[3] += fac0*Real(0.5)*y*(y-Real(1.0)); + stencil(i,j,k)[4] += fac0*Real(0.5)*y*(y+Real(1.0)); } stencil(i,j,k)[0] += fac0*s0; diff --git a/Src/Extern/HYPRE/AMReX_HypreMLABecLap_K.H b/Src/Extern/HYPRE/AMReX_HypreMLABecLap_K.H index ea38bf5037f..129a6a989a8 100644 --- a/Src/Extern/HYPRE/AMReX_HypreMLABecLap_K.H +++ b/Src/Extern/HYPRE/AMReX_HypreMLABecLap_K.H @@ -22,7 +22,7 @@ void hypmlabeclap_mat (GpuArray& sten, int i, int j, in GpuArray, AMREX_SPACEDIM*2> const& bcmsk, GpuArray, AMREX_SPACEDIM*2> const& bcval, GpuArray, AMREX_SPACEDIM*2> const& bcrhs, - int level) + int level, IntVect const& fixed_pt) { Real bxm = b[0] ? b[0](i ,j ,k ) : Real(1.0); Real bxp = b[0] ? b[0](i+1,j ,k ) : Real(1.0); @@ -223,6 +223,12 @@ void hypmlabeclap_mat (GpuArray& sten, int i, int j, in } #endif + + if (fixed_pt == IntVect(AMREX_D_DECL(i,j,k))) { + for (int n = 1; n < 2*AMREX_SPACEDIM+1; ++n) { + sten[n] = Real(0.0); + } + } } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H index 536d4c82b04..783d17bbac0 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H @@ -86,6 +86,8 @@ public: void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, const MF* crse_bcdata=nullptr) override; + void prepareForFluxes (int amrlev, const MF* crse_bcdata = nullptr) override; + void correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b, BCMode bc_mode, const MF* crse_bcdata=nullptr) final; @@ -133,6 +135,10 @@ public: Vector > m_robin_bcval; +#ifdef AMREX_USE_HYPRE + void setInterpBndryHalfWidth (int w) { m_interpbndry_halfwidth = w; } +#endif + protected: bool m_has_metric_term = false; @@ -202,6 +208,8 @@ private: void computeVolInv () const; mutable Vector > m_volinv; // used by solvability fix + + int m_interpbndry_halfwidth = 2; }; template @@ -472,7 +480,9 @@ MLCellLinOpT::defineBC () bc_data.setVal(0.0); m_bndry_cor[amrlev]->setBndryValues(*m_crse_cor_br[amrlev], 0, bc_data, 0, 0, ncomp, - IntVect(this->m_amr_ref_ratio[amrlev-1])); + IntVect(this->m_amr_ref_ratio[amrlev-1]), + InterpBndryDataT::IBD_max_order_DEF, + m_interpbndry_halfwidth); Vector > bclohi (ncomp,Array{{AMREX_D_DECL(BCType::Dirichlet, @@ -544,7 +554,9 @@ MLCellLinOpT::setLevelBC (int amrlev, const MF* a_levelbcdata, const MF* rob m_crse_sol_br[amrlev]->setVal(RT(0.0)); } m_bndry_sol[amrlev]->setBndryValues(*m_crse_sol_br[amrlev], 0, - bcdata, 0, 0, ncomp, br_ref_ratio); + bcdata, 0, 0, ncomp, br_ref_ratio, + InterpBndryDataT::IBD_max_order_DEF, + m_interpbndry_halfwidth); br_ref_ratio = this->m_coarse_data_crse_ratio; } else @@ -639,7 +651,9 @@ MLCellLinOpT::updateSolBC (int amrlev, const MF& crse_bcdata) const m_crse_sol_br[amrlev]->copyFrom(crse_bcdata, 0, 0, 0, ncomp, this->m_geom[amrlev-1][0].periodicity()); m_bndry_sol[amrlev]->updateBndryValues(*m_crse_sol_br[amrlev], 0, 0, ncomp, - IntVect(this->m_amr_ref_ratio[amrlev-1])); + IntVect(this->m_amr_ref_ratio[amrlev-1]), + InterpBndryDataT::IBD_max_order_DEF, + m_interpbndry_halfwidth); } template @@ -652,7 +666,9 @@ MLCellLinOpT::updateCorBC (int amrlev, const MF& crse_bcdata) const m_crse_cor_br[amrlev]->copyFrom(crse_bcdata, 0, 0, 0, ncomp, this->m_geom[amrlev-1][0].periodicity()); m_bndry_cor[amrlev]->updateBndryValues(*m_crse_cor_br[amrlev], 0, 0, ncomp, - IntVect(this->m_amr_ref_ratio[amrlev-1])); + IntVect(this->m_amr_ref_ratio[amrlev-1]), + InterpBndryDataT::IBD_max_order_DEF, + m_interpbndry_halfwidth); } template @@ -1210,6 +1226,15 @@ MLCellLinOpT::solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, MF::Xpay(resid, RT(-1.0), b, 0, 0, ncomp, IntVect(0)); } +template +void +MLCellLinOpT::prepareForFluxes (int amrlev, const MF* crse_bcdata) +{ + if (crse_bcdata != nullptr) { + updateSolBC(amrlev, *crse_bcdata); + } +} + template void MLCellLinOpT::correctionResidual (int amrlev, int mglev, MF& resid, MF& x, const MF& b, diff --git a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H index 03da0874e79..3cc623b761e 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLLinOp.H @@ -376,6 +376,8 @@ public: virtual void solutionResidual (int amrlev, MF& resid, MF& x, const MF& b, const MF* crse_bcdata=nullptr) = 0; + virtual void prepareForFluxes (int /*amrlev*/, const MF* /*crse_bcdata*/ = nullptr) {} + /** * \brief Compute residual for the residual-correction form, resid = b - L(x) * diff --git a/Src/LinearSolvers/MLMG/AMReX_MLMG.H b/Src/LinearSolvers/MLMG/AMReX_MLMG.H index 5078383d5c0..78b2ffdd3df 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLMG.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLMG.H @@ -164,6 +164,8 @@ public: void setHypreStrongThreshold (Real t) noexcept {hypre_strong_threshold = t;} #endif + void prepareForFluxes (Vector const& a_sol); + template void prepareForSolve (Vector const& a_sol, Vector const& a_rhs); @@ -538,6 +540,16 @@ MLMGT::solve (const Vector& a_sol, const Vector& a_rhs, return composite_norminf; } +template +void +MLMGT::prepareForFluxes (Vector const& a_sol) +{ + for (int alev = finest_amr_lev; alev >= 0; --alev) { + const MF* crse_bcdata = (alev > 0) ? a_sol[alev-1] : nullptr; + linop.prepareForFluxes(alev, crse_bcdata); + } +} + template template void