Allow outflow to be zero when inflow is zero while enforcing solvabil…

…ity (#159)

This fixes a bug in the existing code which requires the outflow to be
non-zero even if the inflow is zero for the direction-dependent
boundaries. It should be okay for the outflow to be zero when inflow is
zero, e.g. for finer levels that are not touching the boundary at all,
or even in general.

EBGodunov/hydro_ebgodunov_bcs_K.H

@@ -66,7 +66,7 @@ void SetXBCs (const int i, const int j, const int k, const int n,
             //
             // Note that this is only relevant for 3D, as only 3D ever needs to
             // set the BC beyond the domain face.
-            // This is potentially tricky because the code re-uses some of the
+            // This is potentially tricky because the code reuses some of the
             // space holding the holding the states passed in here; it puts the
             // upwinded intermediate edge state in Imx (which supplies hi)
             // With GPU, I think it's undertermined which Imx(i+..) you'll get here,

Utils/hydro_enforce_inout_solvability.cpp

@@ -6,8 +6,6 @@
 
 using namespace amrex;
 
-constexpr amrex::Real eps = 1.0e-6;
-
 namespace HydroUtils {
 
 namespace {
@@ -320,12 +318,15 @@ void enforceInOutSolvability (
         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);
-        AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
-            outflux > eps,
-            "Cannot enforce solvability, no outflow from the direction dependent boundaries");
 
-        const Real alpha_fcf = influx/outflux;  // flux correction factor
-        correct_outflow(lev, vels_vec, inout_masks, bc_type, domain, alpha_fcf, include_bndry_corners);
+        if ((influx > small_vel) && (outflux < small_vel)) {
+            Abort("Cannot enforce solvability, no outflow from the direction dependent boundaries");
+        } else if ((influx < small_vel) && (outflux < small_vel)) {
+            return; // do nothing
+        } else {
+            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
 }