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

amr-wind/equation_systems/icns/icns_advection.cpp

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

amr-wind/equation_systems/vof/volume_fractions.H

@@ -64,6 +64,66 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal(
     mz = mm1 - mm2;
 }
 
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void youngs_fd_normal_neumann(
+    int i,
+    int j,
+    int k,
+    int ibdy,
+    int jbdy,
+    int kbdy,
+    amrex::Array4<amrex::Real const> const& volfrac,
+    amrex::Real& mx,
+    amrex::Real& my,
+    amrex::Real& mz) noexcept
+{
+    amrex::Real mm1, mm2;
+
+    // Do neumann condition via indices
+    int im1 = ibdy == -1 ? i : i - 1;
+    int jm1 = jbdy == -1 ? j : j - 1;
+    int km1 = kbdy == -1 ? k : k - 1;
+    int ip1 = ibdy == +1 ? i : i + 1;
+    int jp1 = jbdy == +1 ? j : j + 1;
+    int kp1 = kbdy == +1 ? k : k + 1;
+
+    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);
+    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,

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();
+    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_hs_pgrad{true};
+
+    // Coupling options
+    bool m_disable_nodal_proj{false};
+    bool m_disable_mac_proj{false};
+    bool m_replace_gp{false};
+
+    // Verbosity
+    int m_verbose{0};
+
+    // Reinitialization parameters
+    int m_niterations{10};
+    int m_calcconvint{1}; // calctolniter, cconv
+    amrex::Real m_tol = 1e-12;
+    amrex::Real m_rlscale = 1.5;
+    amrex::Real m_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;
 };