AMReX-Fluids · WeiqunZhang · Mar 11, 2024 · Mar 4, 2024 · Mar 4, 2024 · Mar 11, 2024
diff --git a/Docs/source/EB/EB_example.fig b/Docs/source/EB/EB_example.fig
diff --git a/Docs/source/EB/EB_example.png b/Docs/source/EB/EB_example.png
diff --git a/Docs/source/EB/areas_and_volumes.fig b/Docs/source/EB/areas_and_volumes.fig
diff --git a/Docs/source/EB/areas_and_volumes.png b/Docs/source/EB/areas_and_volumes.png
diff --git a/Docs/source/EB/eb_fluxes.fig b/Docs/source/EB/eb_fluxes.fig
diff --git a/Docs/source/EB/eb_fluxes.png b/Docs/source/EB/eb_fluxes.png
diff --git a/Docs/source/EB/redist.fig b/Docs/source/EB/redist.fig
diff --git a/Docs/source/EB/redist.png b/Docs/source/EB/redist.png
diff --git a/Docs/source/Godunov.rst b/Docs/source/Godunov.rst
@@ -231,7 +231,7 @@ Godunov with Embedded Boundaries (EBGodunov)
 
 AMReX-Hydro contains an embedded boundary (EB) aware version of the Godunov algorithm
 discussed above, although with fewer options available.
-This EB implementation employs a piecewise linear method (PLM) :cite:`needref`, and
+This EB implementation employs a piecewise linear method (PLM), and
 always includes any forcing terms *after* the computation of the transverse terms.
 EBGodunov attempts to use fourth-order limited slopes wherever possible, as described in :ref:`EBslopes`.
 

diff --git a/Docs/source/Introduction.rst b/Docs/source/Introduction.rst
@@ -29,7 +29,7 @@ is as follows:
 5. Difference fluxes to create the convective term.
 
 AMReX-Hydro provides routines to support all of these steps.
-Here we group the AMReX-Hydro routines into a few general categories and map them to the step(s) they address:
+Here we group the routines into a few general categories and map them to the step(s) they address:
 
 * :ref:`Schemes`: the fundamental algorithm is either a Method-of-Lines (MOL) or Godunov approach
   (used in steps 1 and 3).
@@ -43,7 +43,7 @@ Here we group the AMReX-Hydro routines into a few general categories and map the
     (not generally used as part of computing the convective term, but used in application codes to define a
     cell-centered velocity update that approximately obeys a divergence constraint).
 
-* :ref:`Redistribution:Redistribution` schemes: to address the "small cell problem" associated with explicit cut
+* :ref:`Redistribution <amrex:sec:EB:redistribution>`: AMReX contains schemes to address the "small cell problem" associated with explicit cut
   cell algorithms (part of step 4 for problems with embedded boundaries).
 
 * :ref:`utilities`: to do things like compute slopes, create fluxes from face-centered values, and

diff --git a/Docs/source/MOL.rst b/Docs/source/MOL.rst
@@ -14,7 +14,7 @@ These alogrithms are applied in the MOL namespace. For API documentation, see
 
 
 Pre-MAC (API ref. `MOL::ExtrapVelToFaces <https://amrex-fluids.github.io/amrex-hydro/Doxygen/html/namespaceMOL.html#acdde2acf756048b8ef0bca332e4bf748>`_)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 For computing the pre-MAC edge states to be MAC-projected, we define on
 every x-face:
@@ -45,7 +45,7 @@ We similarly compute :math:`v_{i,j-\frac{1}{2},k}` on y-faces and
 
 
 Post-MAC (API ref. `MOL::ComputeEdgeState <https://amrex-fluids.github.io/amrex-hydro/Doxygen/html/namespaceMOL.html#acdde2acf756048b8ef0bca332e4bf748>`_)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 Once we have the MAC-projected velocities, we extrapolate all quantities to
 faces as above:
@@ -83,7 +83,7 @@ All slope computations use second-order limited slopes as described in :ref:`EBs
 
 
 Pre-MAC (API ref. `EBMOL::ExtrapVelToFaces <https://amrex-fluids.github.io/amrex-hydro/Doxygen/html/namespaceEBMOL.html#a7add53a153ade9c5cb83e79a61ad1929>`_)
-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 For computing the pre-MAC edge states to be MAC-projected, we define on every x-face with non-zero area fraction:
 

diff --git a/Docs/source/Utilities.rst b/Docs/source/Utilities.rst
@@ -34,8 +34,6 @@ Note that EB is only an option for Cartesian geometries as of this writing.
 
 .. include::   Slopes.rst
 
-.. include::   bcs.rst
-
 .. include::   Fluxes.rst
 
 .. include::   AdvectiveTerm.rst

diff --git a/Docs/source/bcs.rst b/Docs/source/bcs.rst
@@ -1,16 +1,55 @@
+.. role:: cpp(code)
+   :language: c++
+
 .. _bcs:
 
 
 Boundary conditions
--------------------
+====================
 
 AMReX-Hydro uses underlying AMReX functionality in implementing boundary conditions
 (see AMReX's documentation section :ref:`amrex:sec:basics:boundary`).
 Physical boundary conditions, such as
-inflow, outflow, slip/no-slip walls, etc., and are ultimately linked to
+inflow, outflow, slip/no-slip walls, etc. are ultimately linked to
 mathematical Dirichlet or Neumann conditions.
 See ``amrex/Src/Base/AMReX_BC_TYPES.H`` for common physical and mathematical types.
 
+
+.. _mixedBC:
+
+Mixed Boundary Conditions
+--------------------------
+
+.. note::
+
+   This has only been tested in geometries where the Neumann and Dirchlet areas are separated
+   by an embedded boundary, and is not guaranteed to work for other cases.
+
+In addition to the bc types listed in AMReX, mixed boundary conditions can be implemented.
+To this end, the (EB)Godunov advection routines support position-dependent boundary conditions.
+The routines accept a boundary condition MultiFab (BC MF), or Array4. The BC MF is a cell-centered
+(integer) :cpp:`iMultiFab` that carries an :cpp:`amrex::BCType` in the first ghost cell
+and must fully specify the BC on all faces.
+If a position-dependent BCs are passed in, they take precedent and single BC per face :cpp:`BCRecs` are
+ignored.
+
+
+The MacProjector supports mixed boundary conditions by making use of the underlying linear solver's
+Robin BC (:ref:`amrex:sec:linearsolver:bc`) option. Robin boundary conditions are formulated as
+:math:`a\phi + b\frac{\partial\phi}{\partial n} = f`.
+:math:`a`, :math:`b`, and :math:`f` are each a cell-centered (real) :cpp:`MultiFab` that carries
+the relevant values in the first ghost cell. All three must be specified on each level with a call to
+:cpp:`setLevelBC` (see :ref:`amrex:sec:linearsolver:bc` for usage).
+
+The NodalProjector can support mixed boundary conditions through the use of an overset mask
+(see :ref:`amrex:sec:linearsolver:bc`).
+The overset mask specifies a Dirichlet BC with 0 (meaning no solve is needed since the solution is known) or Neumann with 1 (meaning do the solve). Note this is an integer (not bool) MultiFab, so the values must be only either 0 or 1.
+
+
+
+Advective BC details
+--------------------
+
 Domain boundary conditions affect the pre-MAC extrapolated velocities in three ways.
 
 #. Potential impact to the slope computation in cells

diff --git a/Docs/source/conf.py b/Docs/source/conf.py
@@ -107,7 +107,6 @@ def get_amrex_hydro_version():
                     'MOL.rst',
                     'AdvectiveTerm.rst',
                     'Fluxes.rst',
-                    'bcs.rst',
                     'Slopes.rst'
 ]
 

diff --git a/Docs/source/index.rst b/Docs/source/index.rst
@@ -23,6 +23,7 @@ applications.
    Installation Guide <Install>
    Advection Schemes <Schemes>
    Projections
+   Boundary Conditions <bcs>
    Helper functions <Utilities>
 
 .. toctree::

diff --git a/Projections/hydro_MacProjector.H b/Projections/hydro_MacProjector.H
@@ -90,6 +90,8 @@ public:
 
     //! Update Bcoeffs for the linear operator
     void updateBeta (const amrex::Vector<amrex::Array<amrex::MultiFab const*,AMREX_SPACEDIM> >&);
+    //! Update Bcoeff and if Robin BC also reset scalars in order to reuse the linear operator
+    void updateCoeffs (const amrex::Vector<amrex::Array<amrex::MultiFab const*,AMREX_SPACEDIM> >&);
 
 #ifndef AMREX_USE_EB
     void initProjector (amrex::Vector<amrex::BoxArray> const& a_grids,
@@ -184,6 +186,7 @@ private:
 
     bool m_needs_domain_bcs = true;
     amrex::Vector<int> m_needs_level_bcs;
+    bool m_has_robin = false;
 
     // Location of umac -- face center vs face centroid
     amrex::MLMG::Location m_umac_loc;

diff --git a/Projections/hydro_MacProjector.cpp b/Projections/hydro_MacProjector.cpp
@@ -147,15 +147,39 @@ void MacProjector::updateBeta (
         m_poisson == nullptr,
         "MacProjector::updateBeta: should not be called for constant beta");
 
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+        m_has_robin == false,
+        "MacProjector::updateBeta: should not be called with Robin BC. Call updateCoeffs");
+
+    updateCoeffs(a_beta);
+}
+
+void MacProjector::updateCoeffs (
+    const Vector<Array<MultiFab const*, AMREX_SPACEDIM>>& a_beta)
+{
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+        m_linop != nullptr,
+        "MacProjector::updateCoeffs: initProjector must be called before calling this method");
+
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+        m_poisson == nullptr,
+        "MacProjector::updateCoeffs: should not be called for constant beta");
+
     const auto nlevs = int(a_beta.size());
 #ifdef AMREX_USE_EB
     const bool has_eb = a_beta[0][0]->hasEBFabFactory();
     if (has_eb) {
+        if (m_has_robin) {
+            m_eb_abeclap->setScalars(0.0, 1.0);
+        }
         for (int ilev=0; ilev < nlevs; ++ilev)
             m_eb_abeclap->setBCoeffs(ilev, a_beta[ilev], m_beta_loc);
     } else
 #endif
     {
+        if (m_has_robin) {
+            m_abeclap->setScalars(0.0, 1.0);
+        }
         for (int ilev=0; ilev < nlevs; ++ilev)
             m_abeclap->setBCoeffs(ilev, a_beta[ilev]);
     }
@@ -212,9 +236,7 @@ MacProjector::setLevelBC (int amrlev, const MultiFab* levelbcdata, const MultiFa
                                      "setDomainBC must be called before setLevelBC");
     m_linop->setLevelBC(amrlev, levelbcdata, robin_a, robin_b, robin_f);
     m_needs_level_bcs[amrlev] = false;
-    if (robin_a) {
-        m_needs_init = true; // Solver is not safe for reuse with Robin BC
-    }
+    if (robin_a) { m_has_robin = true; }
 }