Multiphase Hybrid - FY23 Q4 milestone (reinitialization/sharpening and pressure approach)

amr-wind/equation_systems/icns/icns_advection.cpp

@@ -55,6 +55,10 @@ MacProjOp::MacProjOp(
 {
     amrex::ParmParse pp("incflo");
     pp.query("density", m_rho_0);
+    amrex::ParmParse pp_ovst("Overset");
+    bool disable_ovst_mac = false;
+    pp_ovst.query("disable_coupled_mac_proj", disable_ovst_mac);
+    m_has_overset = m_has_overset && !disable_ovst_mac;
 }
 
 void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept

amr-wind/equation_systems/vof/volume_fractions.H

@@ -16,27 +16,27 @@ namespace amr_wind::multiphase {
  *  to zero.
  */
 
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
-    int i,
-    int j,
-    int k,
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_finite_difference_normal(
+    const int i,
+    const int j,
+    const int k,
     amrex::Array4<amrex::Real const> const& volfrac,
     amrex::Real& mx,
     amrex::Real& my,
     amrex::Real& mz) noexcept
 {
-    amrex::Real mm1, mm2;
-
-    mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
-          volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
-          2.0 * (volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
-                 volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1)) +
-          4.0 * volfrac(i - 1, j, k);
-    mm2 = volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
-          volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
-          2.0 * (volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
-                 volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1)) +
-          4.0 * volfrac(i + 1, j, k);
+    amrex::Real mm1 =
+        volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
+        volfrac(i - 1, j + 1, k - 1) + volfrac(i - 1, j + 1, k + 1) +
+        2.0 * (volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
+               volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1)) +
+        4.0 * volfrac(i - 1, j, k);
+    amrex::Real mm2 =
+        volfrac(i + 1, j - 1, k - 1) + volfrac(i + 1, j - 1, k + 1) +
+        volfrac(i + 1, j + 1, k - 1) + volfrac(i + 1, j + 1, k + 1) +
+        2.0 * (volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
+               volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1)) +
+        4.0 * volfrac(i + 1, j, k);
     mx = mm1 - mm2;
 
     mm1 = volfrac(i - 1, j - 1, k - 1) + volfrac(i - 1, j - 1, k + 1) +
@@ -64,25 +64,82 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
     mz = mm1 - mm2;
 }
 
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void
+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<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);
+    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);
+    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);
+    mz = mm1 - mm2;
+}
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
-    int i,
-    int j,
-    int k,
+    const int i,
+    const int j,
+    const int k,
     amrex::Array4<amrex::Real const> const& volfrac,
     amrex::Real& mx,
     amrex::Real& my,
     amrex::Real& mz) noexcept
 {
     amrex::Array2D<amrex::Real, 0, AMREX_SPACEDIM + 1, 0, AMREX_SPACEDIM> m;
-    amrex::Real m1, m2;
     // write the plane as: sgn(mx) X =  my Y +  mz Z + alpha
     //                           m00 X = m01 Y + m02 Z + alpha
-    m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
-         volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
-         volfrac(i - 1, j, k);
-    m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
-         volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
-         volfrac(i + 1, j, k);
+    amrex::Real m1 = volfrac(i - 1, j, k - 1) + volfrac(i - 1, j, k + 1) +
+                     volfrac(i - 1, j - 1, k) + volfrac(i - 1, j + 1, k) +
+                     volfrac(i - 1, j, k);
+    amrex::Real m2 = volfrac(i + 1, j, k - 1) + volfrac(i + 1, j, k + 1) +
+                     volfrac(i + 1, j - 1, k) + volfrac(i + 1, j + 1, k) +
+                     volfrac(i + 1, j, k);
 
     m(0, 0) = (m1 > m2) ? 1.0 : -1.0;
 
@@ -180,7 +237,8 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
     }
 
     // Youngs-CIAM scheme */
-    youngs_fd_normal(i, j, k, volfrac, m(3, 0), m(3, 1), m(3, 2));
+    youngs_finite_difference_normal(
+        i, j, k, volfrac, m(3, 0), m(3, 1), m(3, 2));
     // normalize the set (mx,my,mz): |mx|+|my|+|mz| = 1
     constexpr amrex::Real tiny = 1e-20;
     t0 = std::abs(m(3, 0)) + std::abs(m(3, 1)) + std::abs(m(3, 2)) + tiny;
@@ -214,7 +272,10 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void mixed_youngs_central_normal(
  * given that m1+m2+m3=1 (m1,m2,m3>0) and the volumetric fraction volF
  */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(
-    amrex::Real b1, amrex::Real b2, amrex::Real b3, amrex::Real volF) noexcept
+    const amrex::Real b1,
+    const amrex::Real b2,
+    const amrex::Real b3,
+    const amrex::Real volF) noexcept
 {
     using namespace amrex;
 
@@ -296,15 +357,15 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real volume_intercept(
  * (5) calculate volume (NOTE: it is assumed:s0=t0=0; ds0=dt0=1.)
  */
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real cut_volume(
-    amrex::Real m1,
-    amrex::Real m2,
-    amrex::Real m3,
-    amrex::Real alpha,
-    amrex::Real r0,
-    amrex::Real dr0) noexcept
+    const amrex::Real m1,
+    const amrex::Real m2,
+    const amrex::Real m3,
+    const amrex::Real alpha,
+    const amrex::Real r0,
+    const amrex::Real dr0) noexcept
 {
 
-    amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
+    const amrex::Real const_tiny = std::numeric_limits<amrex::Real>::epsilon();
     amrex::Real al;
 
     // move origin to x0
@@ -428,11 +489,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void fit_plane(
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
-    int i,
-    int j,
-    int k,
+    const int i,
+    const int j,
+    const int k,
     amrex::Array4<amrex::Real const> const& volfrac,
-    int n_band = 1) noexcept
+    const int n_band = 1) noexcept
 {
     // n_band must be <= number of vof ghost cells (3)
     constexpr amrex::Real tiny = 1e-12;
@@ -456,14 +517,14 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
 }
 
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real levelset_to_vof(
-    int i,
-    int j,
-    int k,
-    amrex::Real eps,
+    const int i,
+    const int j,
+    const int k,
+    const amrex::Real eps,
     amrex::Array4<amrex::Real const> const& phi) noexcept
 {
     amrex::Real mx, my, mz;
-    youngs_fd_normal(i, j, k, phi, mx, my, mz);
+    youngs_finite_difference_normal(i, j, k, phi, mx, my, mz);
     mx = std::abs(mx / 32.);
     my = std::abs(my / 32.);
     mz = std::abs(mz / 32.);

amr-wind/incflo.cpp

@@ -319,6 +319,9 @@ void incflo::do_advance()
     } else {
         advance();
     }
+    if (m_sim.has_overset()) {
+        m_ovst_ops.post_advance_work();
+    }
 }
 
 // Make a new level from scratch using provided BoxArray and

amr-wind/overset/OversetOps.H

@@ -13,10 +13,52 @@ public:
 
     void initialize(CFDSim& sim);
     void pre_advance_work();
+    void post_advance_work();
 
     void update_gradp();
 
 private:
+    // Functions called within public functions
+    void parameter_output() const;
+    void sharpen_nalu_data();
+    void set_hydrostatic_gradp();
+    void replace_masked_gradp();
+
+    // Check for multiphase sim
+    bool m_vof_exists{false};
+    // Check for perturbational pressure
+    // (will be removed soon)
+    bool m_perturb_p{false};
+
+    // This is the only option for now
+    const bool m_use_hydrostatic_gradp{true};
+
+    // Coupling options
+    bool m_disable_nodal_proj{false};
+    bool m_disable_mac_proj{false};
+    bool m_replace_gradp_postsolve{false};
+
+    // Verbosity
+    int m_verbose{0};
+
+    // Reinitialization parameters
+    int m_n_iterations{10};
+    int m_calc_convg_interval{1}; // calctolniter, cconv
+    amrex::Real m_convg_tol = 1e-12;
+    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
+
+    // Tolerance for VOF-related bound checks
+    const amrex::Real m_vof_tol = 1e-12;
+    // Small number for approximate signed distance function
+    const amrex::Real m_asdf_tiny = 1e-12;
+
+    // Pointer for pressure gradient copy field
+    amr_wind::Field* m_gp_copy{nullptr};
+    // Pointer for MultiPhase physics
+    amr_wind::MultiPhase* m_mphase{nullptr};
+
     CFDSim* m_sim_ptr;
 };