Enforcing inflow-outflow boundary solvability (#129)

Introduces routines to calculate outflux and influx at the in-out
boundaries and then correct the outflux to match with influx.

Three steps:
1. set masks for tagging cells at the boundary for inflow or outflow
2. calculate influx and outflux using reductions
3. correct outflow velocities using the correction factor

---------

Co-authored-by: Ann Almgren <[email protected]>

.gitignore

@@ -54,3 +54,5 @@ Tutorials/EB/Donut/Exec/ls_plt*
 
 # local config
 Tools/GNUMake/Make.local
+
+.vscode*

Docs/source/InOutSolvability.rst

@@ -0,0 +1,28 @@
+.. include:: CustomCommands.rst
+
+Enforcing inflow-outflow solvability
+------------------------------------
+
+This routine enforces solvability for inflow-outflow boundaries,
+which have both inflow and outflow cells.
+A flux correction factor, :math:`\alpha_\text{fcf}` is introduced
+which scales the outflow velocities to match the inflow:
+
+.. math::
+
+    \sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area} + \alpha_\text{fcf} \sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area} = 0
+
+
+The new flux-conserving velocities to be used for the MAC/nodal projections,
+:math:`{\bf u}_\text{fc}`, are calculated from the correction factor as:
+
+.. math::
+    \alpha_\text{fcf} = \frac{-\sum_{\partial\Omega_\text{in}} {\bf u} \cdot {\bf \area}}{\sum_{\partial\Omega_\text{out}} {\bf u} \cdot {\bf \area}},
+
+.. math::
+    {\bf u}_\text{fc} = \alpha_\text{fcf} \cdot {\bf u}, \ \forall {\bf x} \in  \partial\Omega_\text{out}.
+
+It must be noted that this routine currently only accounts for boundaries
+with the math BC ``BCType::direction_dependent``, which is to be used for
+an inflow-outflow boundary. It does not compute or correct the boundary velocities
+over other math BC types such as those representing pure inflow or pure outflow.

Docs/source/Utilities.rst

@@ -38,3 +38,4 @@ Note that EB is only an option for Cartesian geometries as of this writing.
 
 .. include::   AdvectiveTerm.rst
 
+.. include::   InOutSolvability.rst

Utils/CMakeLists.txt

@@ -14,4 +14,5 @@ target_sources(
    hydro_utils.cpp
    hydro_constants.H
    hydro_bcs_K.H
+   hydro_enforce_inout_solvability.cpp
    )

Utils/Make.package

@@ -1,6 +1,7 @@
 CEXE_sources += hydro_utils.cpp
 CEXE_sources += hydro_compute_edgestate_and_flux.cpp
 CEXE_sources += hydro_extrap_vel_to_faces.cpp
+CEXE_sources += hydro_enforce_inout_solvability.cpp
 CEXE_headers += hydro_bcs_K.H
 CEXE_headers += hydro_utils.H
 

Utils/hydro_enforce_inout_solvability.cpp

@@ -0,0 +1,300 @@
+/** \addtogroup Utilities
+ * @{
+ */
+
+#include <hydro_utils.H>
+
+using namespace amrex;
+
+namespace HydroUtils {
+
+namespace {
+
+void set_inout_masks(
+    const int lev,
+    const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
+    Array<iMultiFab, AMREX_SPACEDIM>& inout_masks,
+    const BCRec* bc_type,
+    const Box& domain,
+    const bool corners)
+{
+    for (OrientationIter oit; oit != nullptr; ++oit) {
+        const auto ori = oit();
+        const auto side = ori.faceDir();
+        const int dir = ori.coordDir();
+        const auto islow = ori.isLow();
+        const auto ishigh = ori.isHigh();
+
+        // Multifab for normal velocity
+        const auto& vel_mf = vels_vec[lev][dir];
+
+        // mask iMF for the respective velocity direction
+        auto& inout_mask = inout_masks[dir];
+
+        IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir);
+        // domain extent indices for the velocities
+        int dlo;
+        if (dir_index_type == IndexType::CellIndex::CELL) {
+            // lower boundary is at -1 for cell-centered velocity
+            dlo = domain.smallEnd(dir) - 1;
+        } else {
+            // lower boundary is at  0 for face-centered velocity
+            dlo = domain.smallEnd(dir);
+        }
+        int dhi = domain.bigEnd(dir) + 1;
+
+        // get BCs for the normal velocity and set the boundary index
+        // based on low or high side
+        const BCRec ibcrec = bc_type[dir];
+        int bc, bndry;
+        if (side == Orientation::low) {
+            bc = ibcrec.lo(dir);
+            bndry = dlo;
+        } else {
+            bc = ibcrec.hi(dir);
+            bndry = dhi;
+        }
+
+        // limit influx/outflux calculations to the in-out boundaries only
+        if (bc == BCType::direction_dependent) {
+            for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+                Box box = mfi.validbox();
+
+                // include ghost cells along normal velocity for cell-centered
+                // not for face-centered as boundary lies in valid region
+                if (dir_index_type == IndexType::CellIndex::CELL) {
+                    box.grow(dir, 1);
+                }
+
+                // include boundary corners if specified
+                if (corners) {
+                    box.grow((dir+1)%AMREX_SPACEDIM, 1);
+                    box.grow((dir+2)%AMREX_SPACEDIM, 1);
+                }
+
+                // Enter further only if the box bndry is at the domain bndry
+                if ((islow && (box.smallEnd(dir) == dlo))
+                 || (ishigh && (box.bigEnd(dir) == dhi))) {
+
+                    // create a 2D box normal to dir at the low/high bndry
+                    Box box2d(box); box2d.setRange(dir, bndry);
+
+                    auto vel_arr = vel_mf->array(mfi);
+                    auto inout_mask_arr = inout_mask.array(mfi);
+
+                    // tag cells as inflow or outflow by checking vel direction
+                    ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                    {
+                        if ((side == Orientation::low && vel_arr(i,j,k) >= 0)
+                         || (side == Orientation::high && vel_arr(i,j,k) <= 0)) {
+                            inout_mask_arr(i,j,k) = -1;
+                        } else {
+                            inout_mask_arr(i,j,k) = +1;
+                        }
+                    });
+                }
+            }
+        }
+    }
+}
+
+void compute_influx_outflux(
+    const int lev,
+    const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
+    const Array<iMultiFab, AMREX_SPACEDIM>& inout_masks,
+    const Real* a_dx,
+    Real& influx,
+    Real& outflux,
+    const bool corners)
+{
+    influx = 0.0, outflux = 0.0;
+
+    for (int idim = 0; idim < AMREX_SPACEDIM; idim++) {
+
+        // normal face area
+        const Real ds =
+            a_dx[(idim+1) % AMREX_SPACEDIM] * a_dx[(idim+2) % AMREX_SPACEDIM];
+
+        // Multifab for normal velocity
+        const auto& vel_mf = vels_vec[lev][idim];
+
+        // grow in the respective direction if vel is cell-centered
+        IndexType index_type = vel_mf->ixType();
+        index_type.flip(idim); IntVect ngrow = index_type.ixType();
+
+        // grow in the transverse direction to include boundary corners
+        if (corners) {
+            ngrow[(idim+1)%AMREX_SPACEDIM] = 1;
+            ngrow[(idim+2)%AMREX_SPACEDIM] = 1;
+        }
+
+        // mask iMF for the respective velocity direction
+        const auto& inout_mask = inout_masks[idim];
+
+        // define "multi-arrays" and perform reduction using the mask
+        auto const& vel_ma = vel_mf->const_arrays();
+        auto const& inout_mask_ma = inout_mask.const_arrays();
+
+        influx += ds *
+            ParReduce(TypeList<ReduceOpSum>{},
+                     TypeList<Real>{},
+                     *vel_mf, ngrow,
+           [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
+               noexcept -> GpuTuple<Real>
+           {
+              if (inout_mask_ma[box_no](i,j,k) == -1) {
+                   return { std::abs(vel_ma[box_no](i,j,k)) };
+               } else {
+                   return { 0. };
+               }
+           });
+
+        outflux += ds *
+            ParReduce(TypeList<ReduceOpSum>{},
+                     TypeList<Real>{},
+                     *vel_mf, ngrow,
+           [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
+               noexcept -> GpuTuple<Real>
+           {
+               if (inout_mask_ma[box_no](i,j,k) == 1) {
+                   return { std::abs(vel_ma[box_no](i,j,k)) };
+               } else {
+                   return { 0. };
+               }
+           });
+    }
+    ParallelDescriptor::ReduceRealSum(influx);
+    ParallelDescriptor::ReduceRealSum(outflux);
+}
+
+void correct_outflow(
+    const int lev,
+    const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
+    const Array<iMultiFab, AMREX_SPACEDIM>& inout_masks,
+    const BCRec* bc_type,
+    const Box& domain,
+    const Real alpha_fcf,
+    const bool corners)
+{
+    for (OrientationIter oit; oit != nullptr; ++oit) {
+        const auto ori = oit();
+        const auto side = ori.faceDir();
+        const int dir = ori.coordDir();
+        const auto islow = ori.isLow();
+        const auto ishigh = ori.isHigh();
+
+        // Multifab for normal velocity
+        const auto& vel_mf = vels_vec[lev][dir];
+
+        // mask iMF for the respective velocity direction
+        const auto& inout_mask = inout_masks[dir];
+
+        IndexType::CellIndex dir_index_type = (vel_mf->ixType()).ixType(dir);
+        // domain extent indices for the velocities
+        int dlo;
+        if (dir_index_type == IndexType::CellIndex::CELL) {
+            dlo = domain.smallEnd(dir) - 1; // cell-centered boundary
+        } else {
+            dlo = domain.smallEnd(dir);     // face-centered boundary
+        }
+        int dhi = domain.bigEnd(dir) + 1;
+
+        // get BCs for the normal velocity and set the boundary index
+        const BCRec ibcrec = bc_type[dir];
+        int bc, bndry;
+        if (side == Orientation::low) {
+            bc = ibcrec.lo(dir);
+            bndry = dlo;
+        } else {
+            bc = ibcrec.hi(dir);
+            bndry = dhi;
+        }
+
+        if (bc == BCType::direction_dependent) {
+            for (MFIter mfi(*vel_mf, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+
+                Box box = mfi.validbox();
+                if (dir_index_type == IndexType::CellIndex::CELL) {
+                    box.grow(dir, 1);
+                }
+                if (corners) {
+                    box.grow((dir+1)%AMREX_SPACEDIM, 1);
+                    box.grow((dir+2)%AMREX_SPACEDIM, 1);
+                }
+
+                // Enter further only if the box boundary is at the domain boundary
+                if ((islow && (box.smallEnd(dir) == dlo))
+                 || (ishigh && (box.bigEnd(dir) == dhi))) {
+
+                    // create a 2D box normal to dir at the low/high boundary
+                    Box box2d(box); box2d.setRange(dir, bndry);
+
+                    auto vel_arr = vel_mf->array(mfi);
+                    auto inout_mask_arr = inout_mask.array(mfi);
+
+                    ParallelFor(box2d, [=] AMREX_GPU_DEVICE (int i, int j, int k)
+                    {
+                        if (inout_mask_arr(i,j,k) == 1) {
+                            vel_arr(i,j,k) *= alpha_fcf;
+                        }
+                    });
+                }
+            }
+        }
+    }
+}
+
+} // file-local namespace
+
+void enforceInOutSolvability (
+    const Vector<Array<MultiFab*, AMREX_SPACEDIM>>& vels_vec,
+    const BCRec* bc_type,
+    const Vector<Geometry>& geom,
+    bool include_bndry_corners
+)
+{
+    const Box domain = geom[0].Domain();
+
+    const auto nlevs = int(vels_vec.size());
+    for (int lev = 0; lev < nlevs; ++lev) {
+
+        // masks to tag inflow/outflow at the boundaries
+        // separate iMultifab for each velocity direction
+        //  0 for interior cells,
+        // -1 for inflow bndry cells, +1 for outflow bndry cells
+        Array<iMultiFab, AMREX_SPACEDIM> inout_masks;
+
+        // defining the mask iMultifabs in each direction
+        for (int idim = 0; idim < AMREX_SPACEDIM; idim++)
+        {
+            const auto& vel_mf = vels_vec[lev][idim];    // normal velocity multifab
+
+            // grow in the respective direction if vel is cell-centered
+            // to include the boundary cells
+            IndexType index_type = vel_mf->ixType();
+            index_type.flip(idim); IntVect ngrow = index_type.ixType();
+
+            // grow in the transverse direction to include boundary corners
+            if (include_bndry_corners) {
+                ngrow[(idim+1)%AMREX_SPACEDIM] = 1;
+                ngrow[(idim+2)%AMREX_SPACEDIM] = 1;
+            }
+
+            inout_masks[idim].define(vel_mf->boxArray(), vel_mf->DistributionMap(), 1, ngrow);
+            inout_masks[idim].setVal(0);
+        }
+
+        set_inout_masks(lev, vels_vec, inout_masks, bc_type, domain, include_bndry_corners);
+
+        const Real* a_dx = geom[lev].CellSize();
+        Real influx = 0.0, outflux = 0.0;
+        compute_influx_outflux(lev, vels_vec, inout_masks, a_dx, influx, outflux, include_bndry_corners);
+
+        const Real alpha_fcf = influx/outflux;  // flux correction factor
+        correct_outflow(lev, vels_vec, inout_masks, bc_type, domain, alpha_fcf, include_bndry_corners);
+
+    }   // levels loop
+}
+
+}

Utils/hydro_utils.H

@@ -20,6 +20,17 @@
  */
 
 namespace HydroUtils {
+
+/**
+ * \brief Enforces solvablity by scaling outflow to match with inflow.
+ *
+ */
+void enforceInOutSolvability (
+    const amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>>& vels_vec,
+    amrex::BCRec const* bc_type,
+    const amrex::Vector<amrex::Geometry>& geom,
+    bool include_bndry_corners = false);
+
 /**
  * \brief Compute edge state and flux. Most general version for use with multilevel synchronization.
  *        All other versions ultimately call this one.