Enhance sharpening process for two-phase overset

amr-wind/equation_systems/icns/icns_advection.H

@@ -61,6 +61,7 @@ private:
 #endif
     MLMGOptions m_options;
     bool m_has_overset{false};
+    bool m_disable_mphase_overset_ops{false};
     bool m_need_init{true};
     bool m_variable_density{false};
     bool m_mesh_mapping{false};

amr-wind/equation_systems/icns/icns_advection.cpp

@@ -62,6 +62,11 @@ MacProjOp::MacProjOp(
         pp.query("use_fft", m_use_fft);
     }
 #endif
+
+    // Always off for single-phase flows, this can turn off multiphase overset
+    // operations for all cases
+    amrex::ParmParse pp_o("Overset");
+    pp_o.query("disable_mphase_ops", m_disable_mphase_overset_ops);
 }
 
 void MacProjOp::enforce_inout_solvability(
@@ -78,7 +83,7 @@ void MacProjOp::enforce_inout_solvability(
 void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept
 {
     // Prepare masking for projection
-    if (m_has_overset) {
+    if (m_has_overset && !m_disable_mphase_overset_ops) {
         amr_wind::overset_ops::prepare_mask_cell_for_mac(m_repo);
     }
 
@@ -107,7 +112,7 @@ void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept
 void MacProjOp::init_projector(const amrex::Real beta) noexcept
 {
     // Prepare masking for projection
-    if (m_has_overset) {
+    if (m_has_overset && !m_disable_mphase_overset_ops) {
         amr_wind::overset_ops::prepare_mask_cell_for_mac(m_repo);
     }
 
@@ -365,7 +370,7 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt)
     }
 
     // Prepare masking for remainder of algorithm
-    if (m_has_overset) {
+    if (m_has_overset && !m_disable_mphase_overset_ops) {
         amr_wind::overset_ops::revert_mask_cell_after_mac(m_repo);
     }
 

amr-wind/equation_systems/vof/volume_fractions.H

@@ -5,6 +5,16 @@
 #include <cmath>
 #include "amr-wind/utilities/constants.H"
 
+namespace {
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int
+idx(const int i_off, const int j_off, const int k_off)
+{
+    return (k_off + 1) + 3 * (j_off + 1) + 9 * (i_off + 1);
+}
+
+} // namespace
+
 namespace amr_wind::multiphase {
 
 /** Young's finite-difference gradient scheme.
@@ -70,56 +80,106 @@ youngs_finite_difference_normal_neumann(
     const int i,
     const int j,
     const int k,
-    const int ibdy,
-    const int jbdy,
-    const int kbdy,
+    amrex::Array4<int const> const& iblank,
     amrex::Array4<amrex::Real const> const& volfrac,
     amrex::Real& mx,
     amrex::Real& my,
     amrex::Real& mz) noexcept
 {
-    // Do neumann condition via indices
-    const int im1 = ibdy == -1 ? i : i - 1;
-    const int jm1 = jbdy == -1 ? j : j - 1;
-    const int km1 = kbdy == -1 ? k : k - 1;
-    const int ip1 = ibdy == +1 ? i : i + 1;
-    const int jp1 = jbdy == +1 ? j : j + 1;
-    const int kp1 = kbdy == +1 ? k : k + 1;
-
-    amrex::Real mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
-                      volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
-                      2.0 * (volfrac(im1, jm1, k) + volfrac(im1, jp1, k) +
-                             volfrac(im1, j, km1) + volfrac(im1, j, kp1)) +
-                      4.0 * volfrac(im1, j, k);
-    amrex::Real mm2 = volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
-                      volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
-                      2.0 * (volfrac(ip1, jm1, k) + volfrac(ip1, jp1, k) +
-                             volfrac(ip1, j, km1) + volfrac(ip1, j, kp1)) +
-                      4.0 * volfrac(ip1, j, k);
+    // Do neumann condition using iblank array
+    // Cells that are not iblank == -1 get calculated with extrapolation using
+    // the average of all valid (relying only on iblank == -1 cells) gradients
+    const int ibl_c = iblank(i, j, k);
+    amrex::GpuArray<amrex::Real, 6> dvf{0., 0., 0., 0., 0., 0.};
+    const amrex::IntVect iv{i, j, k};
+    for (int ddir = 0; ddir < 3; ++ddir) {
+        int ct_0 = 0;
+        int ct_1 = 0;
+        const int odir0 = (ddir + 1) % 3;
+        const int odir1 = (ddir + 2) % 3;
+        for (int n0 = -1; n0 <= 1; ++n0) {
+            for (int n1 = -1; n1 <= 1; ++n1) {
+                amrex::IntVect iv_c = iv;
+                iv_c[odir0] += n0;
+                iv_c[odir1] += n1;
+                amrex::IntVect iv_p = iv_c;
+                amrex::IntVect iv_m = iv_c;
+                iv_p[ddir] += 1;
+                iv_m[ddir] += -1;
+                if (iblank(iv_m) == iblank(iv_c) && iblank(iv_m) == ibl_c) {
+                    dvf[ddir] += volfrac(iv_c) - volfrac(iv_m);
+                    ct_0++;
+                }
+                if (iblank(iv_p) == iblank(iv_c) && iblank(iv_p) == ibl_c) {
+                    dvf[ddir + 3] += volfrac(iv_p) - volfrac(iv_c);
+                    ct_1++;
+                }
+            }
+        }
+        dvf[ddir] /= (amrex::Real)ct_0 + constants::EPS;
+        dvf[ddir + 3] /= (amrex::Real)ct_1 + constants::EPS;
+    }
+
+    amrex::GpuArray<amrex::Real, 27> vf;
+    vf.fill(0.);
+    for (int n0 = -1; n0 <= 1; ++n0) {
+        for (int n1 = -1; n1 <= 1; ++n1) {
+            for (int n2 = -1; n2 <= 1; ++n2) {
+                amrex::IntVect iv_l = iv;
+                iv_l[0] += n0;
+                iv_l[1] += n1;
+                iv_l[2] += n2;
+                vf[idx(n0, n1, n2)] =
+                    iblank(iv_l) == ibl_c
+                        ? volfrac(iv_l)
+                        : std::max(
+                              0.,
+                              std::min(
+                                  1., volfrac(iv) +
+                                          (amrex::Real)n0 *
+                                              dvf[0 + 3 * std::max(0, n0)] +
+                                          (amrex::Real)n1 *
+                                              dvf[1 + 3 * std::max(0, n1)] +
+                                          (amrex::Real)n2 *
+                                              dvf[2 + 3 * std::max(0, n2)]));
+            }
+        }
+    }
+
+    amrex::Real mm1 = vf[idx(-1, -1, -1)] + vf[idx(-1, -1, 1)] +
+                      vf[idx(-1, 1, -1)] + vf[idx(-1, 1, 1)] +
+                      2.0 * (vf[idx(-1, -1, 0)] + vf[idx(-1, 1, 0)] +
+                             vf[idx(-1, 0, -1)] + vf[idx(-1, 0, 1)]) +
+                      4.0 * vf[idx(-1, 0, 0)];
+    amrex::Real mm2 = vf[idx(1, -1, -1)] + vf[idx(1, -1, 1)] +
+                      vf[idx(1, 1, -1)] + vf[idx(1, 1, 1)] +
+                      2.0 * (vf[idx(1, -1, 0)] + vf[idx(1, 1, 0)] +
+                             vf[idx(1, 0, -1)] + vf[idx(1, 0, 1)]) +
+                      4.0 * vf[idx(1, 0, 0)];
     mx = mm1 - mm2;
 
-    mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jm1, kp1) +
-          volfrac(ip1, jm1, km1) + volfrac(ip1, jm1, kp1) +
-          2.0 * (volfrac(im1, jm1, k) + volfrac(ip1, jm1, k) +
-                 volfrac(i, jm1, km1) + volfrac(i, jm1, kp1)) +
-          4.0 * volfrac(i, jm1, k);
-    mm2 = volfrac(im1, jp1, km1) + volfrac(im1, jp1, kp1) +
-          volfrac(ip1, jp1, km1) + volfrac(ip1, jp1, kp1) +
-          2.0 * (volfrac(im1, jp1, k) + volfrac(ip1, jp1, k) +
-                 volfrac(i, jp1, km1) + volfrac(i, jp1, kp1)) +
-          4.0 * volfrac(i, jp1, k);
+    mm1 = vf[idx(-1, -1, -1)] + vf[idx(-1, -1, 1)] + vf[idx(1, -1, -1)] +
+          vf[idx(1, -1, 1)] +
+          2.0 * (vf[idx(-1, -1, 0)] + vf[idx(1, -1, 0)] + vf[idx(0, -1, -1)] +
+                 vf[idx(0, -1, 1)]) +
+          4.0 * vf[idx(0, -1, 0)];
+    mm2 = vf[idx(-1, 1, -1)] + vf[idx(-1, 1, 1)] + vf[idx(1, 1, -1)] +
+          vf[idx(1, 1, 1)] +
+          2.0 * (vf[idx(-1, 1, 0)] + vf[idx(1, 1, 0)] + vf[idx(0, 1, -1)] +
+                 vf[idx(0, 1, 1)]) +
+          4.0 * vf[idx(0, 1, 0)];
     my = mm1 - mm2;
 
-    mm1 = volfrac(im1, jm1, km1) + volfrac(im1, jp1, km1) +
-          volfrac(ip1, jm1, km1) + volfrac(ip1, jp1, km1) +
-          2.0 * (volfrac(im1, j, km1) + volfrac(ip1, j, km1) +
-                 volfrac(i, jm1, km1) + volfrac(i, jp1, km1)) +
-          4.0 * volfrac(i, j, km1);
-    mm2 = volfrac(im1, jm1, kp1) + volfrac(im1, jp1, kp1) +
-          volfrac(ip1, jm1, kp1) + volfrac(ip1, jp1, kp1) +
-          2.0 * (volfrac(im1, j, kp1) + volfrac(ip1, j, kp1) +
-                 volfrac(i, jm1, kp1) + volfrac(i, jp1, kp1)) +
-          4.0 * volfrac(i, j, kp1);
+    mm1 = vf[idx(-1, -1, -1)] + vf[idx(-1, 1, -1)] + vf[idx(1, -1, -1)] +
+          vf[idx(1, 1, -1)] +
+          2.0 * (vf[idx(-1, 0, -1)] + vf[idx(1, 0, -1)] + vf[idx(0, -1, -1)] +
+                 vf[idx(0, 1, -1)]) +
+          4.0 * vf[idx(0, 0, -1)];
+    mm2 = vf[idx(-1, -1, 1)] + vf[idx(-1, 1, 1)] + vf[idx(1, -1, 1)] +
+          vf[idx(1, 1, 1)] +
+          2.0 * (vf[idx(-1, 0, 1)] + vf[idx(1, 0, 1)] + vf[idx(0, -1, 1)] +
+                 vf[idx(0, 1, 1)]) +
+          4.0 * vf[idx(0, 0, 1)];
     mz = mm1 - mm2;
 }
 

amr-wind/overset/OversetOps.H

@@ -27,6 +27,10 @@ private:
     // Check for multiphase sim
     bool m_vof_exists{false};
 
+    // Option to turn off all multiphase-specific treatments (for testing and
+    // comparison)
+    bool m_disable_mphase_ops{false};
+
     // Coupling options
     bool m_replace_gradp_postsolve{true};
 
@@ -40,6 +44,8 @@ private:
     amrex::Real m_relative_length_scale = 1.5;
     amrex::Real m_upw_margin = 0.1;
     amrex::Real m_target_cutoff = 0.0; // proc_tgvof_tol
+    // Maximum pseudoCFL
+    amrex::Real m_pCFL = 1.0;
 
     // Tolerance for VOF-related bound checks
     const amrex::Real m_vof_tol = 1e-12;