Add 2D Cartesian/RZ EB Neumann support to MLNodeLaplacian

Docs/sphinx_documentation/source/EB.rst

@@ -633,7 +633,7 @@ AMReX supports multi-level
 1) cell-centered solvers with homogeneous Neumann, homogeneous Dirichlet,
 or inhomogeneous Dirichlet boundary conditions on the EB faces, and
 2) nodal solvers with homogeneous Neumann boundary conditions,
-or inflow velocity conditions on the EB faces.
+2D inhomogeneous Neumann data, or inflow velocity conditions on the EB faces.
 
 To use a cell-centered solver with EB, one builds a linear operator
 :cpp:`MLEBABecLap` with :cpp:`EBFArrayBoxFactory` (instead of a :cpp:`MLABecLaplacian`)
@@ -653,6 +653,12 @@ The usage of this EB-specific class is essentially the same as
 
 The default boundary condition on EB faces is homogeneous Neumann.
 
+For 2D node-based EB solves, :cpp:`MLNodeLaplacian` can set inhomogeneous
+Neumann data on EB faces with :cpp:`setEBInhomogNeumannFlux`. The supplied
+data is the physical flux :math:`\sigma \partial \phi / \partial n_{EB}` at the
+EB surface. In cylindrical RZ solves, the supplied data should not include the
+radial metric factor; the node-based operator applies that factor internally.
+
 To set homogeneous Dirichlet boundary conditions, call
 
 .. highlight:: c++

Docs/sphinx_documentation/source/LinearSolvers.rst

@@ -490,8 +490,28 @@ AMReX supports multi-level solvers for use with embedded boundaries.
 These include
 1) cell-centered solvers with homogeneous Neumann, homogeneous Dirichlet,
 or inhomogeneous Dirichlet boundary conditions on the EB faces, and
-2) nodal solvers with homogeneous Neumann boundary conditions,
-or inflow velocity conditions on the EB faces.
+2) nodal solvers with homogeneous Neumann boundary conditions, 2D
+inhomogeneous Neumann data, or inflow velocity conditions on the EB faces.
+
+For 2D node-based EB solves, :cpp:`MLNodeLaplacian` can set inhomogeneous
+Neumann data on EB faces with
+
+.. highlight:: c++
+
+::
+
+    ml_node_laplacian.setEBInhomogNeumannFlux(lev, eb_neumann_flux);
+
+The first component of :cpp:`eb_neumann_flux` is a cell-centered cut-cell
+MultiFab containing the physical flux
+:math:`\sigma \partial \phi / \partial n_{EB}` at the EB surface. The AMReX EB
+normal :math:`n_{EB}` points from the valid computational region toward the
+covered region. For cylindrical RZ solves, :cpp:`eb_neumann_flux` should not
+include the radial metric factor; :cpp:`MLNodeLaplacian` applies that factor
+internally when forming the EB contribution to the right-hand side. Physical
+domain boundary conditions are still set separately with :cpp:`setDomainBC` and
+related linear-operator calls. The same EB normal convention is used when
+:cpp:`setEBInflowVelocity` forms :math:`u \cdot n_{EB}` on EB faces.
 
 To use a cell-centered solver with EB, one builds a linear operator
 :cpp:`MLEBABecLap` with :cpp:`EBFArrayBoxFactory` (instead of a :cpp:`MLABecLaplacian`)

Src/LinearSolvers/MLMG/AMReX_MLNodeLap_2D_K.H

@@ -1414,7 +1414,8 @@ void mlndlap_res_cf_contrib_cs (int i, int j, int, Array4<Real> const& res,
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void mlndlap_set_stencil (Box const& bx, Array4<Real> const& sten,
                           Array4<Real const> const& sigma,
-                          GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
+                          GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
+                          bool is_rz = false) noexcept
 {
     Real facx = Real(1.0/6.0)*dxinv[0]*dxinv[0];
     Real facy = Real(1.0/6.0)*dxinv[1]*dxinv[1];
@@ -1427,6 +1428,11 @@ void mlndlap_set_stencil (Box const& bx, Array4<Real> const& sten,
         sten(i,j,k,1) = f2xmy*(sigma(i,j-1,k)+sigma(i,j,k));
         sten(i,j,k,2) = fmx2y*(sigma(i-1,j,k)+sigma(i,j,k));
         sten(i,j,k,3) = fxy*sigma(i,j,k);
+        if (is_rz) {
+            Real fp = facy / static_cast<Real>(2*i+1);
+            Real fm = facy / static_cast<Real>(2*i-1);
+            sten(i,j,k,2) += fm*sigma(i-1,j,k) - fp*sigma(i,j,k);
+        }
     });
 }
 
@@ -1883,7 +1889,8 @@ void mlndlap_set_connection (int i, int j, int, Array4<Real> const& conn,
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void mlndlap_set_stencil_eb (int i, int j, int, Array4<Real> const& sten,
                              Array4<Real const> const& sig, Array4<Real const> const& conn,
-                             GpuArray<Real,AMREX_SPACEDIM> const& dxinv) noexcept
+                             GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
+                             bool is_rz = false) noexcept
 {
     Real facx = Real(1./6.)*dxinv[0]*dxinv[0];
     Real facy = Real(1./6.)*dxinv[1]*dxinv[1];
@@ -1893,6 +1900,12 @@ void mlndlap_set_stencil_eb (int i, int j, int, Array4<Real> const& sten,
     sten(i,j,0,2) = Real(2.)*facy*(sig(i-1,j,0)*conn(i-1,j,0,5)+sig(i,j,0)*conn(i,j,0,3))
                             -facx*(sig(i-1,j,0)*conn(i-1,j,0,1)+sig(i,j,0)*conn(i,j,0,1));
     sten(i,j,0,3) = (facx*conn(i,j,0,1)+facy*conn(i,j,0,4))*sig(i,j,0);
+    if (is_rz) {
+        Real fp = facy / static_cast<Real>(2*i+1);
+        Real fm = facy / static_cast<Real>(2*i-1);
+        sten(i,j,0,2) += fm*sig(i-1,j,0)*conn(i-1,j,0,5)
+                       - fp*sig(i  ,j,0)*conn(i  ,j,0,3);
+    }
 }
 
 
@@ -1985,6 +1998,70 @@ void add_eb_flow_contribution (int i, int j, int, Array4<Real> const& rhs,
     }
 }
 
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+Real eb_neumann_surface_term (int i, int j, Real weight,
+                              Array4<Real const> const& bcent,
+                              Array4<Real const> const& eb_neumann_flux,
+                              GpuArray<Real,AMREX_SPACEDIM> const& dx,
+                              GpuArray<Real,AMREX_SPACEDIM> const& problo,
+                              bool is_rz) noexcept
+{
+    if (weight == Real(0.0)) {
+        return Real(0.0);
+    }
+
+    Real metric = Real(1.0);
+    if (is_rz) {
+        metric = problo[0] + (static_cast<Real>(i) + Real(0.5) + bcent(i,j,0,0))*dx[0];
+    }
+
+    return metric * eb_neumann_flux(i,j,0) * weight;
+}
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+void add_eb_neumann_contribution (int i, int j, int, Array4<Real> const& rhs,
+                      Array4<int const> const& msk, GpuArray<Real,AMREX_SPACEDIM> const& dxinv,
+                      GpuArray<Real,AMREX_SPACEDIM> const& dx,
+                      GpuArray<Real,AMREX_SPACEDIM> const& problo,
+                      Array4<Real const> const& bareaarr,
+                      Array4<Real const> const& bcent,
+                      Array4<Real const> const& sintg,
+                      Array4<Real const> const& eb_neumann_flux,
+                      bool is_rz) noexcept
+{
+    using namespace nodelap_detail;
+
+    // eb_neumann_flux is the physical sigma*dphi/dn_EB, with n_EB pointing
+    // from valid cells toward covered cells. RZ metric weighting is applied
+    // here when needed.
+    Real fac_eb = Real(0.25) * dxinv[0];
+
+    if (!msk(i,j,0)) {
+        Real const w_im_jm =          bareaarr(i-1,j-1,0)
+                            + Real(2.)*sintg(i-1,j-1,0,i_B_x )
+                            + Real(2.)*sintg(i-1,j-1,0,i_B_y )
+                            + Real(4.)*sintg(i-1,j-1,0,i_B_xy);
+        Real const w_i_jm  =          bareaarr(i  ,j-1,0)
+                            - Real(2.)*sintg(i  ,j-1,0,i_B_x )
+                            + Real(2.)*sintg(i  ,j-1,0,i_B_y )
+                            - Real(4.)*sintg(i  ,j-1,0,i_B_xy);
+        Real const w_im_j  =          bareaarr(i-1,j  ,0)
+                            + Real(2.)*sintg(i-1,j  ,0,i_B_x )
+                            - Real(2.)*sintg(i-1,j  ,0,i_B_y )
+                            - Real(4.)*sintg(i-1,j  ,0,i_B_xy);
+        Real const w_i_j   =          bareaarr(i  ,j  ,0)
+                            - Real(2.)*sintg(i  ,j  ,0,i_B_x )
+                            - Real(2.)*sintg(i  ,j  ,0,i_B_y )
+                            + Real(4.)*sintg(i  ,j  ,0,i_B_xy);
+
+        rhs(i,j,0) -= fac_eb*(
+           eb_neumann_surface_term(i-1,j-1,w_im_jm,bcent,eb_neumann_flux,dx,problo,is_rz)
+          +eb_neumann_surface_term(i  ,j-1,w_i_jm ,bcent,eb_neumann_flux,dx,problo,is_rz)
+          +eb_neumann_surface_term(i-1,j  ,w_im_j ,bcent,eb_neumann_flux,dx,problo,is_rz)
+          +eb_neumann_surface_term(i  ,j  ,w_i_j  ,bcent,eb_neumann_flux,dx,problo,is_rz));
+    }
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void mlndlap_mknewu_eb (int i, int j, int, Array4<Real> const& u, Array4<Real const> const& p,
                         Array4<Real const> const& sig, Array4<Real const> const& vfrac,

Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.H

@@ -126,6 +126,7 @@ public :
     void getFluxes (const Vector<MultiFab*>& a_flux,
                             const Vector<MultiFab*>& a_sol) const final;
     void unimposeNeumannBC (int amrlev, MultiFab& rhs) const final;
+    void applyInhomogNeumannTerm (int amrlev, MultiFab& rhs) const final;
     Vector<Real> getSolvabilityOffset (int amrlev, int mglev,
                                        MultiFab const& rhs) const override;
     void fixSolvabilityByOffset (int amrlev, int mglev, MultiFab& rhs,
@@ -151,9 +152,26 @@ public :
     void buildIntegral ();
     void buildSurfaceIntegral ();
 
-    // Tells the solver that EB boundaries have inflow specified by "eb_vel"
+    /**
+     * \brief Set velocity data on embedded boundaries.
+     *
+     * The solver uses `eb_vel dot n_EB` on EB faces, where the AMReX EB normal
+     * `n_EB` points from the valid/fluid region toward the covered/body region.
+     */
     void setEBInflowVelocity (int amrlev, const MultiFab& eb_vel);
 
+    /**
+     * \brief Set inhomogeneous Neumann data on embedded boundaries.
+     *
+     * The first component of `eb_neumann_flux` is the cell-centered cut-cell
+     * value of the physical flux `sigma*dphi/dn_EB` at the EB surface. AMReX
+     * EB normals point from the valid/fluid region toward the covered/body
+     * region. For RZ solves, this flux should not include the radial metric
+     * factor; the solver applies that factor internally. This is currently
+     * implemented for 2D EB nodal solves.
+     */
+    void setEBInhomogNeumannFlux (int amrlev, const MultiFab& eb_neumann_flux);
+
 #endif
 
 #if defined(AMREX_USE_HYPRE) && (AMREX_SPACEDIM > 1)
@@ -209,6 +227,7 @@ private:
     bool m_surface_integral_built = false;
     bool m_build_surface_integral = false;
     Vector<std::unique_ptr<MultiFab> > m_eb_vel_dot_n;
+    Vector<std::unique_ptr<MultiFab> > m_eb_neumann_flux;
 #endif
 
     bool m_use_gauss_seidel     = true;

Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian.cpp

@@ -86,6 +86,7 @@ MLNodeLaplacian::define (const Vector<Geometry>& a_geom,
     m_integral.resize(m_num_amr_levels);
     m_surface_integral.resize(m_num_amr_levels);
     m_eb_vel_dot_n.resize(m_num_amr_levels);
+    m_eb_neumann_flux.resize(m_num_amr_levels);
     for (int amrlev = 0; amrlev < m_num_amr_levels; ++amrlev)
     {
         m_integral[amrlev] = std::make_unique<MultiFab>(m_grids[amrlev][0],
@@ -1079,6 +1080,42 @@ MLNodeLaplacian::setEBInflowVelocity (int amrlev, const MultiFab& eb_vel)
                                                     *m_factory[amrlev][0]);
     // Turn on flag for building surface integrals
     m_build_surface_integral = true;
+    m_surface_integral_built = false;
+}
+
+void
+MLNodeLaplacian::setEBInhomogNeumannFlux (int amrlev, const MultiFab& eb_neumann_flux)
+{
+#if (AMREX_SPACEDIM != 2)
+    amrex::ignore_unused(amrlev, eb_neumann_flux);
+    amrex::Abort("MLNodeLaplacian::setEBInhomogNeumannFlux is currently implemented for 2D only");
+#else
+    const int mglev = 0;
+    const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][mglev].get());
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(factory != nullptr,
+        "MLNodeLaplacian::setEBInhomogNeumannFlux requires an EB factory");
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(eb_neumann_flux.nComp() >= 1,
+        "MLNodeLaplacian::setEBInhomogNeumannFlux requires at least one component");
+
+    if (m_eb_neumann_flux[amrlev] == nullptr || m_eb_neumann_flux[amrlev]->nComp() != 1) {
+        m_eb_neumann_flux[amrlev] = std::make_unique<MultiFab>(
+            m_grids[amrlev][mglev], m_dmap[amrlev][mglev],
+            1, 1, MFInfo(), *m_factory[amrlev][mglev]);
+    }
+
+    m_eb_neumann_flux[amrlev]->setVal(0.0);
+    m_eb_neumann_flux[amrlev]->ParallelCopy(eb_neumann_flux, 0, 0, 1,
+                                            m_geom[amrlev][mglev].periodicity());
+    m_eb_neumann_flux[amrlev]->FillBoundary(m_geom[amrlev][mglev].periodicity());
+
+    const int ncomp_si = 3;
+    m_surface_integral[amrlev] = std::make_unique<MultiFab>(m_grids[amrlev][0],
+                                                    m_dmap[amrlev][0],
+                                                    ncomp_si, 1, MFInfo(),
+                                                    *m_factory[amrlev][0]);
+    m_build_surface_integral = true;
+    m_surface_integral_built = false;
+#endif
 }
 #endif
 

Src/LinearSolvers/MLMG/AMReX_MLNodeLaplacian_misc.cpp

@@ -1059,6 +1059,71 @@ MLNodeLaplacian::compDivergence (const Vector<MultiFab*>& rhs, const Vector<Mult
     compRHS(rhs, vel, Vector<const MultiFab*>(), Vector<MultiFab*>());
 }
 
+void
+MLNodeLaplacian::applyInhomogNeumannTerm (int amrlev, MultiFab& rhs) const
+{
+#if !defined(AMREX_USE_EB)
+    amrex::ignore_unused(amrlev, rhs);
+#elif (AMREX_SPACEDIM != 2)
+    amrex::ignore_unused(amrlev, rhs);
+    if (!m_eb_neumann_flux.empty() && m_eb_neumann_flux[amrlev]) {
+        amrex::Abort("MLNodeLaplacian::applyInhomogNeumannTerm is currently implemented for 2D only");
+    }
+#else
+    if (m_eb_neumann_flux.empty() || m_eb_neumann_flux[amrlev] == nullptr) { return; }
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_surface_integral_built,
+        "MLNodeLaplacian::applyInhomogNeumannTerm requires built EB surface integrals");
+
+    const int mglev = 0;
+    const Geometry& geom = m_geom[amrlev][mglev];
+    const auto dxinvarr = geom.InvCellSizeArray();
+    const auto dxarr = geom.CellSizeArray();
+    const auto probloarr = geom.ProbLoArray();
+    const bool is_rz = m_is_rz;
+
+    const auto *factory = dynamic_cast<EBFArrayBoxFactory const*>(m_factory[amrlev][mglev].get());
+    if (!factory || factory->isAllRegular()) { return; }
+
+    const auto& flags = factory->getMultiEBCellFlagFab();
+    const MultiCutFab& barea = factory->getBndryArea();
+    const MultiCutFab& bcent = factory->getBndryCent();
+    const MultiFab& sintg = *m_surface_integral[amrlev];
+    const MultiFab& eb_neumann_flux = *m_eb_neumann_flux[amrlev];
+    const iMultiFab& dmsk = *m_dirichlet_mask[amrlev][mglev];
+
+    MFItInfo mfi_info;
+    if (Gpu::notInLaunchRegion()) { mfi_info.EnableTiling().SetDynamic(true); }
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+    for (MFIter mfi(rhs, mfi_info); mfi.isValid(); ++mfi)
+    {
+        const Box& bx = mfi.tilebox();
+        const auto& flag = flags[mfi];
+        FabType typ = flag.getType(amrex::enclosedCells(bx));
+        if (typ == FabType::singlevalued)
+        {
+            Array4<Real> const& rhsarr = rhs.array(mfi);
+            Array4<int const> const& dmskarr = dmsk.const_array(mfi);
+            Array4<Real const> const& bareaarr = barea.const_array(mfi);
+            Array4<Real const> const& bcarr = bcent.const_array(mfi);
+            Array4<Real const> const& sintgarr = sintg.const_array(mfi);
+            Array4<Real const> const& ebneuarr = eb_neumann_flux.const_array(mfi);
+
+            AMREX_HOST_DEVICE_FOR_3D(bx, i, j, k,
+            {
+                add_eb_neumann_contribution(i, j, k, rhsarr, dmskarr, dxinvarr,
+                                            dxarr, probloarr, bareaarr, bcarr,
+                                            sintgarr, ebneuarr, is_rz);
+            });
+        }
+    }
+
+    nodalSync(amrlev, mglev, rhs);
+#endif
+}
+
 void
 MLNodeLaplacian::compRHS (const Vector<MultiFab*>& rhs, const Vector<MultiFab*>& vel,  // NOLINT(readability-convert-member-functions-to-static)
                           const Vector<const MultiFab*>& rhnd,