From bc88b60eca477582c313fc67e26bf5b6838915ee Mon Sep 17 00:00:00 2001 From: Weiqun Zhang Date: Thu, 23 May 2024 14:13:31 -0700 Subject: [PATCH] InterpBndryData: Make changes for multi-level hypre We want to limit the max width of the interpolation stencil at the coarse/fine interface. --- Src/Boundary/AMReX_InterpBndryData.H | 37 ++++++++++++++-------- Src/Boundary/AMReX_InterpBndryData_1D_K.H | 3 +- Src/Boundary/AMReX_InterpBndryData_2D_K.H | 16 ++++++---- Src/Boundary/AMReX_InterpBndryData_3D_K.H | 6 ++-- Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H | 20 ++++++------ 5 files changed, 49 insertions(+), 33 deletions(-) diff --git a/Src/Boundary/AMReX_InterpBndryData.H b/Src/Boundary/AMReX_InterpBndryData.H index 6dabc8ac481..0a561fae71b 100644 --- a/Src/Boundary/AMReX_InterpBndryData.H +++ b/Src/Boundary/AMReX_InterpBndryData.H @@ -91,7 +91,7 @@ public: */ void setBndryValues (BndryRegisterT const& crse, int c_start, const MF& fine, int f_start, int bnd_start, int num_comp, const IntVect& ratio, - int max_order = IBD_max_order_DEF); + int max_order = IBD_max_order_DEF, int max_width = 2); /** * \brief Update boundary values at coarse/fine boundaries @@ -127,14 +127,15 @@ template void InterpBndryDataT::setPhysBndryValues (const MF& mf, int mf_start, int bnd_start, int num_comp) { - AMREX_ASSERT(this->grids == mf.boxArray()); + AMREX_ASSERT(mf.empty() || this->grids == mf.boxArray()); const Box& fine_domain = this->geom.Domain(); #ifdef AMREX_USE_OMP #pragma omp parallel if (Gpu::notInLaunchRegion()) #endif - for (MFIter mfi(mf,MFItInfo().SetDynamic(true)); mfi.isValid(); ++mfi) + for (MFIter mfi(this->grids, this->DistributionMap(), + MFItInfo().SetDynamic(true)); mfi.isValid(); ++mfi) { const Box& bx = mfi.validbox(); for (OrientationIter fi; fi; ++fi) { @@ -143,14 +144,22 @@ InterpBndryDataT::setPhysBndryValues (const MF& mf, int mf_start, int bnd_st { // Physical bndry, copy from grid. auto & bnd_fab = this->bndry[face][mfi]; - auto const& src_fab = mf[mfi]; auto const& bnd_array = bnd_fab.array(); - auto const& src_array = src_fab.const_array(); - const Box& b = src_fab.box() & bnd_fab.box(); - AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( b, num_comp, i, j, k, n, - { - bnd_array(i,j,k,n+bnd_start) = src_array(i,j,k,n+mf_start); - }); + if (mf.empty()) { + const Box& b = bnd_fab.box(); + AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( b, num_comp, i, j, k, n, + { + bnd_array(i,j,k,n+bnd_start) = value_type(0); + }); + } else { + auto const& src_fab = mf[mfi]; + auto const& src_array = src_fab.const_array(); + const Box& b = src_fab.box() & bnd_fab.box(); + AMREX_HOST_DEVICE_PARALLEL_FOR_4D ( b, num_comp, i, j, k, n, + { + bnd_array(i,j,k,n+bnd_start) = src_array(i,j,k,n+mf_start); + }); + } } } } @@ -160,7 +169,7 @@ template void InterpBndryDataT::setBndryValues (BndryRegisterT const& crse, int c_start, const MF& fine, int f_start, int bnd_start, int num_comp, const IntVect& ratio, - int max_order) + int max_order, int max_width) { AMREX_ASSERT(this->grids == fine.boxArray()); @@ -205,7 +214,7 @@ InterpBndryDataT::setBndryValues (BndryRegisterT const& crse, int c_star { interpbndrydata_x_o3(it,jt,kt,n,bdry_array,bnd_start, crse_array,c_start,rr, - mask_array, is_not_covered); + mask_array, is_not_covered, max_width); }); break; } @@ -216,7 +225,7 @@ InterpBndryDataT::setBndryValues (BndryRegisterT const& crse, int c_star { interpbndrydata_y_o3(it,jt,kt,n,bdry_array,bnd_start, crse_array,c_start,rr, - mask_array, is_not_covered); + mask_array, is_not_covered, max_width); }); break; } @@ -227,7 +236,7 @@ InterpBndryDataT::setBndryValues (BndryRegisterT const& crse, int c_star { interpbndrydata_z_o3(it,jt,kt,n,bdry_array,bnd_start, crse_array,c_start,rr, - mask_array, is_not_covered); + mask_array, is_not_covered, max_width); }); break; } diff --git a/Src/Boundary/AMReX_InterpBndryData_1D_K.H b/Src/Boundary/AMReX_InterpBndryData_1D_K.H index 1a17a7bc21e..6c8e2e56cbc 100644 --- a/Src/Boundary/AMReX_InterpBndryData_1D_K.H +++ b/Src/Boundary/AMReX_InterpBndryData_1D_K.H @@ -22,7 +22,8 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interpbndrydata_x_o3 (int i, int /*j*/, int /*k*/, int n, Array4 const& bdry, int nb, Array4 const& crse, int nc, Dim3 const& r, - Array4 const& /*mask*/, int /*not_covered*/) noexcept + Array4 const& /*mask*/, int /*not_covered*/, + int /*max_width*/) noexcept { int ic = amrex::coarsen(i,r.x); bdry(i,0,0,n+nb) = crse(ic,0,0,n+nc); diff --git a/Src/Boundary/AMReX_InterpBndryData_2D_K.H b/Src/Boundary/AMReX_InterpBndryData_2D_K.H index f33d81c7639..584636eee95 100644 --- a/Src/Boundary/AMReX_InterpBndryData_2D_K.H +++ b/Src/Boundary/AMReX_InterpBndryData_2D_K.H @@ -23,7 +23,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interpbndrydata_x_o3 (int i, int j, int /*k*/, int n, Array4 const& bdry, int nb, Array4 const& crse, int nc, Dim3 const& r, - Array4 const& mask, int not_covered) noexcept + Array4 const& mask, int not_covered, int max_width) noexcept { int ic = amrex::coarsen(i,r.x); int jc = amrex::coarsen(j,r.y); @@ -37,7 +37,8 @@ void interpbndrydata_x_o3 (int i, int j, int /*k*/, int n, x[NN] = T(-1.0); y[NN] = crse(ic,jc-1,0,n+nc); ++NN; - } else if (mask.contains(i,(jc+2)*r.y,0) && + } else if (max_width >= 2 && + mask.contains(i,(jc+2)*r.y,0) && mask (i,(jc+2)*r.y,0) == not_covered && crse.contains(ic,jc+2,0)) { x[NN] = T(2.0); y[NN] = crse(ic,jc+2,0,n+nc); @@ -48,7 +49,8 @@ void interpbndrydata_x_o3 (int i, int j, int /*k*/, int n, x[NN] = T(1.0); y[NN] = crse(ic,jc+1,0,n+nc); ++NN; - } else if (mask.contains(i,jc*r.y-r.y-1,0) && + } else if (max_width >= 2 && + mask.contains(i,jc*r.y-r.y-1,0) && mask (i,jc*r.y-r.y-1,0) == not_covered && crse.contains(ic,jc-2,0)) { x[NN] = T(-2.0); y[NN] = crse(ic,jc-2,0,n+nc); @@ -73,7 +75,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interpbndrydata_y_o3 (int i, int j, int /*k*/, int n, Array4 const& bdry, int nb, Array4 const& crse, int nc, Dim3 const& r, - Array4 const& mask, int not_covered) noexcept + Array4 const& mask, int not_covered, int max_width) noexcept { int ic = amrex::coarsen(i,r.x); int jc = amrex::coarsen(j,r.y); @@ -87,7 +89,8 @@ void interpbndrydata_y_o3 (int i, int j, int /*k*/, int n, x[NN] = T(-1.0); y[NN] = crse(ic-1,jc,0,n+nc); ++NN; - } else if (mask.contains((ic+2)*r.x,j,0) && + } else if (max_width >= 2 && + mask.contains((ic+2)*r.x,j,0) && mask ((ic+2)*r.x,j,0) == not_covered && crse.contains(ic+2,jc,0)) { x[NN] = T(2.0); y[NN] = crse(ic+2,jc,0,n+nc); @@ -98,7 +101,8 @@ void interpbndrydata_y_o3 (int i, int j, int /*k*/, int n, x[NN] = T(1.0); y[NN] = crse(ic+1,jc,0,n+nc); ++NN; - } else if (mask.contains(ic*r.x-r.x-1,j,0) && + } else if (max_width >= 2 && + mask.contains(ic*r.x-r.x-1,j,0) && mask (ic*r.x-r.x-1,j,0) == not_covered && crse.contains(ic-2,jc,0)) { x[NN] = T(-2.0); y[NN] = crse(ic-2,jc,0,n+nc); diff --git a/Src/Boundary/AMReX_InterpBndryData_3D_K.H b/Src/Boundary/AMReX_InterpBndryData_3D_K.H index 992ccfcff89..661b3d94c9d 100644 --- a/Src/Boundary/AMReX_InterpBndryData_3D_K.H +++ b/Src/Boundary/AMReX_InterpBndryData_3D_K.H @@ -23,7 +23,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interpbndrydata_x_o3 (int i, int j, int k, int n, Array4 const& bdry, int nb, Array4 const& crse, int nc, Dim3 const& r, - Array4 const& mask, int not_covered) noexcept + Array4 const& mask, int not_covered, int /*max_width*/) noexcept { int ic = amrex::coarsen(i,r.x); int jc = amrex::coarsen(j,r.y); @@ -56,7 +56,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interpbndrydata_y_o3 (int i, int j, int k, int n, Array4 const& bdry, int nb, Array4 const& crse, int nc, Dim3 const& r, - Array4 const& mask, int not_covered) noexcept + Array4 const& mask, int not_covered, int /*max_width*/) noexcept { int ic = amrex::coarsen(i,r.x); int jc = amrex::coarsen(j,r.y); @@ -90,7 +90,7 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void interpbndrydata_z_o3 (int i, int j, int k, int n, Array4 const& bdry, int nb, Array4 const& crse, int nc, Dim3 const& r, - Array4 const& mask, int not_covered) noexcept + Array4 const& mask, int not_covered, int /*max_width*/) noexcept { int ic = amrex::coarsen(i,r.x); int jc = amrex::coarsen(j,r.y); diff --git a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H index 536d4c82b04..311d04be56e 100644 --- a/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H +++ b/Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H @@ -202,6 +202,8 @@ private: void computeVolInv () const; mutable Vector > m_volinv; // used by solvability fix + + int m_interpbndry_max_width; }; template @@ -472,7 +474,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]), + InterpBndryData::IBD_max_order_DEF, + m_interpbndry_max_width); Vector > bclohi (ncomp,Array{{AMREX_D_DECL(BCType::Dirichlet, @@ -505,15 +509,11 @@ MLCellLinOpT::setLevelBC (int amrlev, const MF* a_levelbcdata, const MF* rob const int ncomp = this->getNComp(); - MF zero; IntVect ng(1); if (this->hasHiddenDimension()) { ng[this->hiddenDirection()] = 0; } - if (a_levelbcdata == nullptr) { - zero.define(this->m_grids[amrlev][0], this->m_dmap[amrlev][0], ncomp, ng); - zero.setVal(RT(0.0)); - } else { - AMREX_ALWAYS_ASSERT(a_levelbcdata->nGrowVect().allGE(ng)); - } + AMREX_ALWAYS_ASSERT(a_levelbcdata == nullptr || a_levelbcdata->nGrowVect().allGE(ng)); + + MF zero; const MF& bcdata = (a_levelbcdata == nullptr) ? zero : *a_levelbcdata; IntVect br_ref_ratio(-1); @@ -544,7 +544,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, + InterpBndryData::IBD_max_order_DEF, + m_interpbndry_max_width); br_ref_ratio = this->m_coarse_data_crse_ratio; } else