MLEBABecLap: Support Robin BC at Domain Boundaries

Src/LinearSolvers/MLMG/AMReX_MLABecLaplacian.H

@@ -190,15 +190,15 @@ public:
  Array<FAB*,AMREX_SPACEDIM> const& flux,
  FAB const& sol, int face_only, int ncomp);
 
-protected:
-
- bool m_needs_update = true;
-
  RT m_a_scalar = std::numeric_limits<RT>::quiet_NaN();
  RT m_b_scalar = std::numeric_limits<RT>::quiet_NaN();
  Vector<Vector<MF> > m_a_coeffs;
  Vector<Vector<Array<MF,AMREX_SPACEDIM> > > m_b_coeffs;
 
+protected:
+
+ bool m_needs_update = true;
+
  Vector<int> m_is_singular;
 
  [[nodiscard]] bool supportRobinBC () const noexcept override { return true; }
@@ -474,29 +474,29 @@ MLABecLaplacianT<MF>::applyMetricTermsCoeffs ()
 // \tilde{alpha}_i = alpha_i + (1-B) beta_{i+1/2} / h^2
 // \tilde{rhs}_i = rhs_i + A beta_{i+1/2} / h^2
 //
-template <typename MF>
-void
-MLABecLaplacianT<MF>::applyRobinBCTermsCoeffs ()
+namespace detail {
+template <typename LP>
+void applyRobinBCTermsCoeffs (LP& linop)
 {
- if (!(this->hasRobinBC())) { return; }
+ using RT = typename LP::RT;
 
- const int ncomp = this->getNComp();
+ const int ncomp = linop.getNComp();
  bool reset_alpha = false;
- if (m_a_scalar == RT(0.0)) {
- m_a_scalar = RT(1.0);
+ if (linop.m_a_scalar == RT(0.0)) {
+ linop.m_a_scalar = RT(1.0);
  reset_alpha = true;
  }
- const RT bovera = m_b_scalar/m_a_scalar;
+ const RT bovera = linop.m_b_scalar/linop.m_a_scalar;
 
- for (int amrlev = 0; amrlev < this->m_num_amr_levels; ++amrlev) {
+ for (int amrlev = 0; amrlev < linop.NAMRLevels(); ++amrlev) {
  const int mglev = 0;
- const Box& domain = this->m_geom[amrlev][mglev].Domain();
- const RT dxi = static_cast<RT>(this->m_geom[amrlev][mglev].InvCellSize(0));
- const RT dyi = static_cast<RT>((AMREX_SPACEDIM >= 2) ? this->m_geom[amrlev][mglev].InvCellSize(1) : Real(1.0));
- const RT dzi = static_cast<RT>((AMREX_SPACEDIM == 3) ? this->m_geom[amrlev][mglev].InvCellSize(2) : Real(1.0));
+ const Box& domain = linop.Geom(amrlev,mglev).Domain();
+ const RT dxi = static_cast<RT>(linop.Geom(amrlev,mglev).InvCellSize(0));
+ const RT dyi = static_cast<RT>((AMREX_SPACEDIM >= 2) ? linop.Geom(amrlev,mglev).InvCellSize(1) : Real(1.0));
+ const RT dzi = static_cast<RT>((AMREX_SPACEDIM == 3) ? linop.Geom(amrlev,mglev).InvCellSize(2) : Real(1.0));
 
  if (reset_alpha) {
- m_a_coeffs[amrlev][mglev].setVal(RT(0.0));
+ linop.m_a_coeffs[amrlev][mglev].setVal(RT(0.0));
  }
 
  MFItInfo mfi_info;
@@ -505,20 +505,20 @@ MLABecLaplacianT<MF>::applyRobinBCTermsCoeffs ()
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (Gpu::notInLaunchRegion())
 #endif
- for (MFIter mfi(m_a_coeffs[amrlev][mglev], mfi_info); mfi.isValid(); ++mfi)
+ for (MFIter mfi(linop.m_a_coeffs[amrlev][mglev], mfi_info); mfi.isValid(); ++mfi)
  {
  const Box& vbx = mfi.validbox();
- auto const& afab = m_a_coeffs[amrlev][mglev].array(mfi);
+ auto const& afab = linop.m_a_coeffs[amrlev][mglev].array(mfi);
  for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
- auto const& bfab = m_b_coeffs[amrlev][mglev][idim].const_array(mfi);
+ auto const& bfab = linop.m_b_coeffs[amrlev][mglev][idim].const_array(mfi);
  const Box& blo = amrex::adjCellLo(vbx,idim);
  const Box& bhi = amrex::adjCellHi(vbx,idim);
  bool outside_domain_lo = !(domain.contains(blo));
  bool outside_domain_hi = !(domain.contains(bhi));
  if ((!outside_domain_lo) && (!outside_domain_hi)) { continue; }
  for (int icomp = 0; icomp < ncomp; ++icomp) {
- auto const& rbc = (*(this->m_robin_bcval[amrlev]))[mfi].const_array(icomp*3);
- if (this->m_lobc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_lo)
+ auto const& rbc = (*(linop.m_robin_bcval[amrlev]))[mfi].const_array(icomp*3);
+ if (linop.m_lobc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_lo)
  {
  if (idim == 0) {
  RT fac = bovera*dxi*dxi;
@@ -546,7 +546,7 @@ MLABecLaplacianT<MF>::applyRobinBCTermsCoeffs ()
  });
  }
  }
- if (this->m_hibc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_hi)
+ if (linop.m_hibc_orig[icomp][idim] == LinOpBCType::Robin && outside_domain_hi)
  {
  if (idim == 0) {
  RT fac = bovera*dxi*dxi;
@@ -579,6 +579,16 @@ MLABecLaplacianT<MF>::applyRobinBCTermsCoeffs ()
  }
  }
 }
+} // namespace detail
+
+template <typename MF>
+void
+MLABecLaplacianT<MF>::applyRobinBCTermsCoeffs ()
+{
+ if (this->hasRobinBC()) {
+ detail::applyRobinBCTermsCoeffs(*this);
+ }
+}
 
 template <typename MF>
 void

Src/LinearSolvers/MLMG/AMReX_MLCellLinOp.H

@@ -131,6 +131,8 @@ public:
  RT location;
  };
 
+ Vector<std::unique_ptr<MF> > m_robin_bcval;
+
 protected:
 
  bool m_has_metric_term = false;
@@ -182,8 +184,6 @@ protected:
  };
  Vector<Vector<std::unique_ptr<BndryCondLoc> > > m_bcondloc;
 
- Vector<std::unique_ptr<MF> > m_robin_bcval;
-
  // used to save interpolation coefficients of the first interior cells
  mutable Vector<Vector<BndryRegisterT<MF>> > m_undrrelxr;
 

Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.H

@@ -114,6 +114,8 @@ public:
  void getEBFluxes (const Vector<MultiFab*>& a_flux,
  const Vector<MultiFab*>& a_sol) const override;
 
+ void applyRobinBCTermsCoeffs ();
+
 #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
  [[nodiscard]] std::unique_ptr<Hypre> makeHypre (Hypre::Interface hypre_interface) const override;
 #endif
@@ -122,6 +124,11 @@ public:
  [[nodiscard]] std::unique_ptr<PETScABecLap> makePETSc () const override;
 #endif
 
+ Real m_a_scalar = std::numeric_limits<Real>::quiet_NaN();
+ Real m_b_scalar = std::numeric_limits<Real>::quiet_NaN();
+ Vector<Vector<MultiFab> > m_a_coeffs;
+ Vector<Vector<Array<MultiFab,AMREX_SPACEDIM> > > m_b_coeffs;
+
 protected:
 
  int m_ncomp = 1;
@@ -131,10 +138,6 @@ protected:
  Location m_beta_loc; // Location of coefficients: face centers or face centroids
  Location m_phi_loc; // Location of solution variable: cell centers or cell centroids
 
- Real m_a_scalar = std::numeric_limits<Real>::quiet_NaN();
- Real m_b_scalar = std::numeric_limits<Real>::quiet_NaN();
- Vector<Vector<MultiFab> > m_a_coeffs;
- Vector<Vector<Array<MultiFab,AMREX_SPACEDIM> > > m_b_coeffs;
  Vector<Vector<iMultiFab> > m_cc_mask;
 
  Vector<std::unique_ptr<MultiFab> > m_eb_phi;
@@ -154,6 +157,8 @@ protected:
  const Vector<MultiFab*>& b_eb);
  void averageDownCoeffs ();
  void averageDownCoeffsToCoarseAmrLevel (int flev);
+
+ [[nodiscard]] bool supportRobinBC () const noexcept override { return true; }
 };
 
 }

Src/LinearSolvers/MLMG/AMReX_MLEBABecLap.cpp

@@ -685,6 +685,8 @@ MLEBABecLap::prepareForSolve ()
 
  MLCellABecLap::prepareForSolve();
 
+ applyRobinBCTermsCoeffs();
+
  averageDownCoeffs();
 
  if (m_eb_phi[0]) {
@@ -1285,6 +1287,14 @@ MLEBABecLap::getEBFluxes (const Vector<MultiFab*>& a_flux, const Vector<MultiFab
  }
 }
 
+void
+MLEBABecLap::applyRobinBCTermsCoeffs ()
+{
+ if (this->hasRobinBC()) {
+ detail::applyRobinBCTermsCoeffs(*this);
+ }
+}
+
 #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
 std::unique_ptr<Hypre>
 MLEBABecLap::makeHypre (Hypre::Interface hypre_interface) const

Src/LinearSolvers/MLMG/AMReX_MLLinOp.H

@@ -488,6 +488,22 @@ public:
 
  [[nodiscard]] bool isMFIterSafe (int amrlev, int mglev1, int mglev2) const;
 
+ //! Return the number of AMR levels
+ [[nodiscard]] int NAMRLevels () const noexcept { return m_num_amr_levels; }
+
+ //! Return the number of MG levels at given AMR level
+ [[nodiscard]] int NMGLevels (int amrlev) const noexcept { return m_num_mg_levels[amrlev]; }
+
+ [[nodiscard]] const Geometry& Geom (int amr_lev, int mglev=0) const noexcept { return m_geom[amr_lev][mglev]; }
+
+ // BC
+ Vector<Array<BCType, AMREX_SPACEDIM> > m_lobc;
+ Vector<Array<BCType, AMREX_SPACEDIM> > m_hibc;
+ // Need to save the original copy because we change the BC type to
+ // Neumann for inhomogeneous Neumann and Robin.
+ Vector<Array<BCType, AMREX_SPACEDIM> > m_lobc_orig;
+ Vector<Array<BCType, AMREX_SPACEDIM> > m_hibc_orig;
+
 protected:
 
  static constexpr int mg_coarsen_ratio = 2;
@@ -544,14 +560,6 @@ protected:
  };
  std::unique_ptr<CommContainer> m_raii_comm;
 
- // BC
- Vector<Array<BCType, AMREX_SPACEDIM> > m_lobc;
- Vector<Array<BCType, AMREX_SPACEDIM> > m_hibc;
- // Need to save the original copy because we change the BC type to
- // Neumann for inhomogeneous Neumann and Robin.
- Vector<Array<BCType, AMREX_SPACEDIM> > m_lobc_orig;
- Vector<Array<BCType, AMREX_SPACEDIM> > m_hibc_orig;
-
  Array<Real, AMREX_SPACEDIM> m_domain_bloc_lo {{AMREX_D_DECL(0._rt,0._rt,0._rt)}};
  Array<Real, AMREX_SPACEDIM> m_domain_bloc_hi {{AMREX_D_DECL(0._rt,0._rt,0._rt)}};
 
@@ -561,20 +569,12 @@ protected:
  const MF* m_coarse_data_for_bc = nullptr;
  MF m_coarse_data_for_bc_raii;
 
- //! Return the number of AMR levels
- [[nodiscard]] int NAMRLevels () const noexcept { return m_num_amr_levels; }
-
- //! Return the number of MG levels at given AMR level
- [[nodiscard]] int NMGLevels (int amrlev) const noexcept { return m_num_mg_levels[amrlev]; }
-
  //! Return AMR refinement ratios
  [[nodiscard]] const Vector<int>& AMRRefRatio () const noexcept { return m_amr_ref_ratio; }
 
  //! Return AMR refinement ratio at given AMR level
  [[nodiscard]] int AMRRefRatio (int amr_lev) const noexcept { return m_amr_ref_ratio[amr_lev]; }
 
- [[nodiscard]] const Geometry& Geom (int amr_lev, int mglev=0) const noexcept { return m_geom[amr_lev][mglev]; }
-
  [[nodiscard]] FabFactory<FAB> const* Factory (int amr_lev, int mglev=0) const noexcept {
  return m_factory[amr_lev][mglev].get();
  }