Beginnings of multiphase hybrid (Exawind#915)

* option to sharpen nalu vof field or not

* putting previous commit into practice

* important options for sharpening and masking

* able to turn off mac projection

* getting toward a working sharpening
-- sloppy commit

* provides stability for very small epsilon (interface thickness)

* removing unnecessary things, clean up

* more unnecessary

* default to a smaller interface thickness

* masking fluxes is sufficient, no need to mask whole term

* density fluxes and update (to make sure it works)

* best upwinding approach so far

* begin to look at the rest of the iteration

* ability to turn off strong overset coupling to nodal projection

* initialize pressure in sloshing tank case

* modify pressure by hydrostatic source term
-- commenting or removing extra stuff

* momentum sharpening as included in paper
* also, found 1.5 * i_th works for pseudo-CFL

* to go with last change

* ability to specify tolerance for convergence, not just steps

* to go with previous commit

* keep vof bounded

* initialize pressure field for other cases (like waves)

* specify relative length scale by input argument

* move copy to old, misc

* Neumann version of youngs normal
- corresponding version of levelset2vof

* discrete sharpening, much improved from other one

* likely unused stuff, add for history

* more reasonable asdf function, better for small nalu domains

* more flexible flux computation and commented debugging stuff

* populate ghost cells for iblank fields
-huge bug fix for this method

* making vof margin for flux upwinding/downwinding an input argument

* changing defaults of sharpening tool

* commented print statements for debugging

* include boundedness of vof and recalc of density

* new FuzzyInterface case

* calculate face velocity in the upwinding/downwinding form like the vof fluxes

* changing the sharpen margin default
* essential to better results with flat interface

* do density adjustment as intended

* options to ignore nalu pressure and to turn off pressure source term in sharpening

* cleaning up the code and putting in more options

* replacing the pressure gradient with hydrostatic
-- should make this optional

* reapply hydrostatic guess of pressure gradient

* turning latest commits into input options

* switching to full upwinding for velocity sharpening

* ramping the pseudo-time seems to work better
- less noisy velocity field

* ability to sharpen pressure gradient

* 3 big changes
- calculate pseudo dt based on pseudo CFL
- make routine suited to multiple levels
- remove unnecessary parts of  disable_pressure_from_nalu

* set current working version of options as default

amr-wind/core/SimTime.H

@@ -40,6 +40,8 @@ public:
      */
     bool new_timestep();
 
+    void increment_timestep() {++m_time_index;}
+
     /** Return true if simulation should continue
      */
     bool continue_simulation() const;

amr-wind/equation_systems/icns/icns_advection.cpp

@@ -50,6 +50,11 @@ MacProjOp::MacProjOp(
 {
     amrex::ParmParse pp("incflo");
     pp.query("density", m_rho_0);
+    bool disable_omac{true};
+    pp.query("disable_overset_mac", disable_omac);
+    if (m_has_overset && disable_omac) {
+        m_has_overset = false;
+    }
 }
 
 void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept

amr-wind/equation_systems/vof/vof_advection.H

@@ -23,6 +23,9 @@ struct AdvectionOp<VOF, fvm::Godunov>
     {
         amrex::ParmParse pp_multiphase("VOF");
         pp_multiphase.query("remove_debris", m_rm_debris);
+        pp_multiphase.query("sharpen_acq_field", m_sharpen_acq);
+        pp_multiphase.query("sharpen_dest_field", m_sharpen_dest);
+        pp_multiphase.query("replace_masked", m_replace_mask);
 
         // Setup density factor arrays for multiplying velocity flux
         fields_in.repo.declare_face_normal_field(
@@ -34,7 +37,9 @@ struct AdvectionOp<VOF, fvm::Godunov>
         const FieldState /*unused*/,
         const amrex::Real /*unused*/,
         const amrex::Real /*unused*/)
-    {}
+    {
+        // Would be best to sharpen old vof and update old density here
+    }
 
     void operator()(const FieldState /*unused*/, const amrex::Real dt)
     {
@@ -43,13 +48,38 @@ struct AdvectionOp<VOF, fvm::Godunov>
 
         auto& repo = fields.repo;
         const auto& geom = repo.mesh().Geom();
+        const int nlevels = repo.num_active_levels();
 
         auto& aa_x = repo.get_field("advalpha_x");
         auto& aa_y = repo.get_field("advalpha_y");
         auto& aa_z = repo.get_field("advalpha_z");
 
-        // cppcheck-suppress constVariable
-        auto& dof_field = fields.field;
+        // Old and new states
+        auto& dof_old = fields.field.state(amr_wind::FieldState::Old);
+        auto& dof_new = fields.field;
+        // Working state of vof is nph, to keep others untouched during step
+        auto& dof_field = fields.field.state(amr_wind::FieldState::NPH);
+
+        // Sharpen starting vof field if hybrid solver is being used
+        if (repo.int_field_exists("iblank_cell") && m_sharpen_acq) {
+            auto& f_iblank = repo.get_int_field("iblank_cell");
+            amr_wind::multiphase::sharpen_acquired_vof(
+                nlevels, f_iblank, dof_old);
+        }
+        // Sharpen future vof field if hybrid solver is being used
+        if (repo.int_field_exists("iblank_cell") && m_sharpen_dest) {
+            auto& f_iblank = repo.get_int_field("iblank_cell");
+            amr_wind::multiphase::sharpen_acquired_vof(
+                nlevels, f_iblank, dof_new);
+        }
+
+        // Initialize as old state values
+        for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
+            amrex::MultiFab::Copy(
+                dof_field(lev), dof_old(lev), 0, 0, dof_field.num_comp(),
+                dof_field.num_grow());
+        }
+
         //
         // Advect volume using Implicit Eulerian Sweeping method with PLIC
         // reconstruction
@@ -71,8 +101,6 @@ struct AdvectionOp<VOF, fvm::Godunov>
             isweep = 1;
         }
 
-        const int nlevels = repo.num_active_levels();
-
         amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> fluxes(
             nlevels);
         amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> advas(
@@ -86,20 +114,6 @@ struct AdvectionOp<VOF, fvm::Godunov>
             advas[lev][2] = &aa_z(lev);
         }
 
-        // Sharpen acquired vof field if hybrid solver is being used
-        if (repo.int_field_exists("iblank_cell")) {
-            auto& f_iblank = repo.get_int_field("iblank_cell");
-            amr_wind::multiphase::sharpen_acquired_vof(
-                nlevels, f_iblank, dof_field);
-            // Update old vof so that old density can be updated
-            // cppcheck-suppress constVariableReference
-            auto& old_dof_field = dof_field.state(amr_wind::FieldState::Old);
-            for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
-                amrex::MultiFab::Copy(
-                    old_dof_field(lev), dof_field(lev), 0, 0,
-                    dof_field.num_comp(), dof_field.num_grow());
-            }
-        }
         // Split advection step 1, with cmask calculation
         multiphase::split_advection_step(
             isweep, 0, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
@@ -112,6 +126,20 @@ struct AdvectionOp<VOF, fvm::Godunov>
         multiphase::split_advection_step(
             isweep, 2, nlevels, dof_field, fluxes, (*fluxC), advas, u_mac,
             v_mac, w_mac, dof_field.bc_type(), geom, dt, m_rm_debris);
+
+        // Replace masked cells using overset
+        if (repo.int_field_exists("iblank_cell") && m_replace_mask) {
+            auto& f_iblank = repo.get_int_field("iblank_cell");
+            amr_wind::multiphase::replace_masked_vof(
+                nlevels, f_iblank, dof_field, dof_new);
+        }
+
+        // Copy working version of vof to new state
+        for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
+            amrex::MultiFab::Copy(
+                dof_new(lev), dof_field(lev), 0, 0, dof_field.num_comp(),
+                dof_field.num_grow());
+        }
     }
 
     PDEFields& fields;
@@ -120,6 +148,9 @@ struct AdvectionOp<VOF, fvm::Godunov>
     Field& w_mac;
     int isweep = 0;
     bool m_rm_debris{true};
+    bool m_sharpen_acq{false};
+    bool m_sharpen_dest{false};
+    bool m_replace_mask{true};
     // Lagrangian transport is deprecated, only Eulerian is supported
 };
 

amr-wind/equation_systems/vof/vof_hybridsolver_ops.H

@@ -45,6 +45,36 @@ static void sharpen_acquired_vof(
     }
 }
 
+static void replace_masked_vof(
+    const int nlevels,
+    amr_wind::IntField& f_iblank,
+    amr_wind::Field& f_vof,
+    amr_wind::Field& f_vof_new)
+{
+    // Sharpen data from nalu-wind (in iblank regions)
+    for (int lev = 0; lev < nlevels; ++lev) {
+        auto& iblank = f_iblank(lev);
+        auto& vof = f_vof(lev);
+        const auto& vof_new = f_vof_new(lev);
+
+        for (amrex::MFIter mfi(iblank); mfi.isValid(); ++mfi) {
+            const auto& gbx = mfi.growntilebox();
+            const amrex::Array4<const int>& native_flag =
+                iblank.const_array(mfi);
+            const amrex::Array4<amrex::Real>& volfrac = vof.array(mfi);
+            const amrex::Array4<const amrex::Real>& vfmasked =
+                vof_new.const_array(mfi);
+            amrex::ParallelFor(
+                gbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+                    // In iblanked regions, sharpen VOF and limit it
+                    volfrac(i, j, k) = (native_flag(i, j, k) > 0)
+                                           ? volfrac(i, j, k)
+                                           : vfmasked(i, j, k);
+                });
+        }
+    }
+}
+
 } // namespace amr_wind::multiphase
 
 #endif // VOF_HYBRIDSOLVER_OPS.H

amr-wind/equation_systems/vof/volume_fractions.H

@@ -64,6 +64,65 @@ 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.H

@@ -173,10 +173,25 @@ private:
     bool m_do_initial_proj = true;
     int m_initial_iterations = 3;
 
+    // Sharpening multiphase for overset
+    int m_sharpen_iterations = 10;
+    amrex::Real m_sharpen_tolerance = 1e-12;
+    int m_sharpen_calctolniter = 1;
+    amrex::Real m_sharpen_rlscale = 1.5;
+    amrex::Real m_sharpen_margin = 0.1;
+    amrex::Real m_sharpen_proctg_tol = 0.0;
+    bool m_sharpen_hs_pressure = false;
+    bool m_sharpen_hsp_guess = true;
+    bool m_sharpen_hsp_replace = true;
+    bool m_sharpen_gradp = false;
+
     bool m_constant_density = true;
 
     bool m_use_godunov = true;
 
+    // Other choices for overset
+    bool m_disable_onodal = true;
+
     // Prescribe advection velocity
     bool m_prescribe_vel = false;
 

amr-wind/incflo.cpp

@@ -388,6 +388,21 @@ void incflo::init_physics_and_pde()
         // Check for if velocity is prescribed
         amrex::ParmParse pp("incflo");
         pp.query("prescribe_velocity", m_prescribe_vel);
+
+        // Check for sharpening iterations (overset feature)
+        pp.query("sharpen_iterations", m_sharpen_iterations);
+        pp.query("sharpen_tolerance", m_sharpen_tolerance);
+        pp.query("sharpen_calctol_niter", m_sharpen_calctolniter);
+        pp.query("sharpen_rlscale", m_sharpen_rlscale);
+        pp.query("sharpen_margin", m_sharpen_margin);
+        pp.query("sharpen_targetvof_tol", m_sharpen_proctg_tol);
+        pp.query("sharpen_hs_pressure", m_sharpen_hs_pressure);
+        pp.query("sharpen_guess_hsp", m_sharpen_hsp_guess);
+        pp.query("sharpen_replace_hsp", m_sharpen_hsp_replace);
+        pp.query("sharpen_pressure_grad", m_sharpen_gradp);
+
+        // Determine if overset values should be forced into projection
+        pp.query("disable_overset_nodal", m_disable_onodal);
     }
     m_sim.create_turbulence_model();
 

amr-wind/incflo_advance.cpp

@@ -8,6 +8,7 @@
 #include "amr-wind/utilities/console_io.H"
 #include "amr-wind/utilities/PostProcessing.H"
 #include "AMReX_MultiFabUtil.H"
+#include "amr-wind/overset/sharpen_nalu_data.H"
 
 using namespace amrex;
 
@@ -190,11 +191,30 @@ void incflo::ApplyPredictor(bool incremental_projection)
     const auto& density_old = density_new.state(amr_wind::FieldState::Old);
     auto& density_nph = density_new.state(amr_wind::FieldState::NPH);
 
-    // Recalculate incoming pressure gradient field if overset
+    auto gp_copy = density().repo().create_scratch_field(3);
+    auto& gp = density().repo().get_field("gp");
+
+    // Process data for overset multiphase
     if (sim().has_overset()) {
-        UpdateGradP(
-            (density_old).vec_const_ptrs(), m_time.current_time(),
-            m_time.deltaT());
+        // Sharpen nalu fields
+        amr_wind::overset::SharpenNaluDataDiscrete(
+            sim(), m_sharpen_iterations, m_sharpen_tolerance,
+            m_sharpen_calctolniter, m_sharpen_rlscale, m_sharpen_margin,
+            m_sharpen_proctg_tol, m_sharpen_hs_pressure, m_sharpen_gradp);
+        // Recalculate pressure gradient with incoming sharpened field
+        if (!m_sharpen_gradp || sim().time().current_time() == 0.0) {
+            UpdateGradP(
+                (density_old).vec_const_ptrs(), m_time.current_time(),
+                m_time.deltaT());
+        }
+        if (m_sharpen_hsp_guess ||
+            (m_sharpen_gradp && sim().time().current_time() == 0.0)) {
+            amr_wind::overset::ReplaceMaskedGradP(sim());
+        }
+        if (m_sharpen_gradp) {
+            amr_wind::overset::CopyGradP(
+                *gp_copy, gp, sim().repo().num_active_levels());
+        }
     }
 
     // *************************************************************************************
@@ -269,6 +289,10 @@ void incflo::ApplyPredictor(bool incremental_projection)
         }
     }
 
+    if (m_verbose > 2) {
+        PrintMaxVelLocations("before pre advection");
+    }
+
     // Extrapolate and apply MAC projection for advection velocities
     icns().pre_advection_actions(amr_wind::FieldState::Old);
 
@@ -396,6 +420,13 @@ void incflo::ApplyPredictor(bool incremental_projection)
         (density_new).vec_const_ptrs(), new_time, m_time.deltaT(),
         incremental_projection);
 
+    if (m_sharpen_hsp_replace) {
+        amr_wind::overset::ReapplyModifiedGradP(sim());
+    }
+    if (m_sharpen_gradp) {
+        amr_wind::overset::ReapplyOversetGradP(*gp_copy, sim());
+    }
+
     if (m_verbose > 2) {
         PrintMaxVelLocations("after nodal projection");
     }

amr-wind/overset/TiogaInterface.H

@@ -104,6 +104,8 @@ private:
 
     CFDSim& m_sim;
 
+    bool m_disable_pressure_from_nalu{true};
+
     //! IBLANK on cell centered fields
     IntField& m_iblank_cell;