Change to multiphase hybrid coupling - both projections

amr-wind/equation_systems/icns/icns_advection.cpp

@@ -4,6 +4,7 @@
 #include "amr-wind/core/MLMGOptions.H"
 #include "amr-wind/utilities/console_io.H"
 #include "amr-wind/wind_energy/ABL.H"
+#include "amr-wind/overset/overset_ops_routines.H"
 
 #include "AMReX_MultiFabUtil.H"
 #include "hydro_MacProjector.H"
@@ -57,10 +58,6 @@ 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::enforce_inout_solvability(
@@ -76,6 +73,12 @@ void MacProjOp::enforce_inout_solvability(
 
 void MacProjOp::init_projector(const MacProjOp::FaceFabPtrVec& beta) noexcept
 {
+    // Prepare masking for projection
+    if (m_has_overset) {
+        amr_wind::overset_ops::prepare_mask_cell_for_mac(m_repo);
+    }
+
+    // Prepare projector
     m_mac_proj = std::make_unique<Hydro::MacProjector>(
         m_repo.mesh().Geom(0, m_repo.num_active_levels() - 1));
     m_mac_proj->initProjector(
@@ -99,6 +102,12 @@ 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) {
+        amr_wind::overset_ops::prepare_mask_cell_for_mac(m_repo);
+    }
+
+    // Prepare projector
     m_mac_proj = std::make_unique<Hydro::MacProjector>(
         m_repo.mesh().Geom(0, m_repo.num_active_levels() - 1));
     m_mac_proj->initProjector(
@@ -308,6 +317,11 @@ void MacProjOp::operator()(const FieldState fstate, const amrex::Real dt)
         }
     }
 
+    // Prepare masking for remainder of algorithm
+    if (m_has_overset) {
+        amr_wind::overset_ops::revert_mask_cell_after_mac(m_repo);
+    }
+
     io::print_mlmg_info("MAC_projection", m_mac_proj->getMLMG());
 }
 

amr-wind/equation_systems/vof/volume_fractions.H

@@ -3,6 +3,7 @@
 
 #include <AMReX_FArrayBox.H>
 #include <cmath>
+#include "amr-wind/utilities/constants.H"
 
 namespace amr_wind::multiphase {
 
@@ -493,10 +494,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool interface_band(
     const int j,
     const int k,
     amrex::Array4<amrex::Real const> const& volfrac,
-    const int n_band = 1) noexcept
+    const int n_band = 1,
+    const amrex::Real tiny = constants::TIGHT_TOL) noexcept
 {
     // n_band must be <= number of vof ghost cells (3)
-    constexpr amrex::Real tiny = 1e-12;
+
     // Check if near interface
     amrex::Real VOF_max = 0.0;
     amrex::Real VOF_min = 1.0;

amr-wind/overset/OversetOps.H

@@ -22,19 +22,13 @@ private:
     void parameter_output() const;
     void sharpen_nalu_data();
     void form_perturb_pressure();
-    void set_hydrostatic_gradp();
     void replace_masked_gradp();
 
     // Check for multiphase sim
     bool m_vof_exists{false};
 
-    // To avoid using sharpened pressure, use hydrostatic pressure
-    bool m_use_hydrostatic_gradp{false};
-
     // Coupling options
-    bool m_disable_nodal_proj{false};
-    bool m_disable_mac_proj{false};
-    bool m_replace_gradp_postsolve{false};
+    bool m_replace_gradp_postsolve{true};
 
     // Verbosity
     int m_verbose{0};

amr-wind/overset/OversetOps.cpp

@@ -24,12 +24,7 @@ void OversetOps::initialize(CFDSim& sim)
     pp.query("reinit_target_cutoff", m_target_cutoff);
 
     // Queries for coupling options
-    pp.query("use_hydrostatic_gradp", m_use_hydrostatic_gradp);
     pp.query("replace_gradp_postsolve", m_replace_gradp_postsolve);
-    // OversetOps does not control these coupling options, merely reports them
-    pp.query("disable_coupled_nodal_proj", m_disable_nodal_proj);
-    pp.query("disable_coupled_mac_proj", m_disable_mac_proj);
-
     pp.query("verbose", m_verbose);
 
     m_vof_exists = m_sim_ptr->repo().field_exists("vof");
@@ -51,48 +46,48 @@ void OversetOps::initialize(CFDSim& sim)
                    "turned on, but it will remain inactive because "
                    "perturbational pressure is turned off.\n";
         }
-    }
-    if (m_replace_gradp_postsolve) {
-        m_gp_copy = &m_sim_ptr->repo().declare_field("gp_copy", 3);
+        if (m_replace_gradp_postsolve) {
+            m_gp_copy = &m_sim_ptr->repo().declare_field("gp_copy", 3);
+        }
     }
 
     parameter_output();
 }
 
 void OversetOps::pre_advance_work()
 {
-    if (!(m_vof_exists && m_use_hydrostatic_gradp)) {
-        if (m_mphase != nullptr) {
-            // Avoid modifying pressure upon initialization, assume pressure = 0
-            if (m_mphase->perturb_pressure() &&
-                m_sim_ptr->time().current_time() > 0.0) {
-                // Modify to be consistent with internal source terms
-                form_perturb_pressure();
-            }
+    if (m_mphase != nullptr) {
+        // Avoid modifying pressure upon initialization, assume pressure = 0
+        if (m_mphase->perturb_pressure() &&
+            m_sim_ptr->time().current_time() > 0.0) {
+            // Modify to be consistent with internal source terms
+            form_perturb_pressure();
         }
-        // Update pressure gradient using updated overset pressure field
-        update_gradp();
     }
+    // Update pressure gradient using updated overset pressure field
+    update_gradp();
 
     if (m_vof_exists) {
         // Reinitialize fields
         sharpen_nalu_data();
-        if (m_use_hydrostatic_gradp) {
-            // Use hydrostatic pressure gradient
-            set_hydrostatic_gradp();
-        } else {
-            // Update pressure gradient using sharpened pressure field
-            update_gradp();
-        }
-    }
-
-    // If pressure gradient will be replaced, store current pressure gradient
-    if (m_replace_gradp_postsolve) {
-        const auto& gp = m_sim_ptr->repo().get_field("gp");
-        for (int lev = 0; lev < m_sim_ptr->repo().num_active_levels(); ++lev) {
-            amrex::MultiFab::Copy(
-                (*m_gp_copy)(lev), gp(lev), 0, 0, gp(lev).nComp(),
-                (m_gp_copy)->num_grow());
+        // Update pressure gradient using sharpened pressure field
+        update_gradp();
+        // Calculate vof-dependent node mask
+        const auto& iblank = m_sim_ptr->repo().get_int_field("iblank_node");
+        const auto& vof = m_sim_ptr->repo().get_field("vof");
+        auto& mask = m_sim_ptr->repo().get_int_field("mask_node");
+        overset_ops::iblank_node_to_mask_vof(iblank, vof, mask);
+
+        // If pressure gradient will be replaced, store current pressure
+        // gradient
+        if (m_replace_gradp_postsolve) {
+            const auto& gp = m_sim_ptr->repo().get_field("gp");
+            for (int lev = 0; lev < m_sim_ptr->repo().num_active_levels();
+                 ++lev) {
+                amrex::MultiFab::Copy(
+                    (*m_gp_copy)(lev), gp(lev), 0, 0, gp(lev).nComp(),
+                    (m_gp_copy)->num_grow());
+            }
         }
     }
 
@@ -163,7 +158,7 @@ void OversetOps::update_gradp()
 void OversetOps::post_advance_work()
 {
     // Replace and reapply pressure gradient if requested
-    if (m_replace_gradp_postsolve) {
+    if (m_vof_exists && m_replace_gradp_postsolve) {
         replace_masked_gradp();
     }
 }
@@ -175,38 +170,30 @@ void OversetOps::post_advance_work()
 void OversetOps::parameter_output() const
 {
     // Print the details
-    if (m_verbose > 0) {
+    if (m_verbose > 0 && m_vof_exists) {
         // Important parameters
         amrex::Print() << "\nOverset Coupling Parameters: \n"
-                       << "---- Coupled nodal projection : "
-                       << !m_disable_nodal_proj << std::endl
-                       << "---- Coupled MAC projection   : "
-                       << !m_disable_mac_proj << std::endl
                        << "---- Replace overset pres grad: "
                        << m_replace_gradp_postsolve << std::endl;
-        if (m_vof_exists) {
-            amrex::Print() << "---- Perturbational pressure  : "
-                           << m_mphase->perturb_pressure() << std::endl
-                           << "---- Reconstruct true pressure: "
-                           << m_mphase->reconstruct_true_pressure()
-                           << std::endl;
-            amrex::Print() << "Overset Reinitialization Parameters:\n"
-                           << "---- Maximum iterations   : " << m_n_iterations
-                           << std::endl
-                           << "---- Convergence tolerance: " << m_convg_tol
-                           << std::endl
-                           << "---- Relative length scale: "
-                           << m_relative_length_scale << std::endl
-                           << "---- Upwinding VOF margin : " << m_upw_margin
+        amrex::Print() << "---- Perturbational pressure  : "
+                       << m_mphase->perturb_pressure() << std::endl
+                       << "---- Reconstruct true pressure: "
+                       << m_mphase->reconstruct_true_pressure() << std::endl;
+        amrex::Print() << "Overset Reinitialization Parameters:\n"
+                       << "---- Maximum iterations   : " << m_n_iterations
+                       << std::endl
+                       << "---- Convergence tolerance: " << m_convg_tol
+                       << std::endl
+                       << "---- Relative length scale: "
+                       << m_relative_length_scale << std::endl
+                       << "---- Upwinding VOF margin : " << m_upw_margin
+                       << std::endl;
+        if (m_verbose > 1) {
+            // Less important or less used parameters
+            amrex::Print() << "---- Calc. conv. interval : "
+                           << m_calc_convg_interval << std::endl
+                           << "---- Target field cutoff  : " << m_target_cutoff
                            << std::endl;
-            if (m_verbose > 1) {
-                // Less important or less used parameters
-                amrex::Print()
-                    << "---- Calc. conv. interval : " << m_calc_convg_interval
-                    << std::endl
-                    << "---- Target field cutoff  : " << m_target_cutoff
-                    << std::endl;
-            }
         }
         amrex::Print() << std::endl;
     }
@@ -414,35 +401,6 @@ void OversetOps::form_perturb_pressure()
     }
 }
 
-void OversetOps::set_hydrostatic_gradp()
-{
-    const auto& repo = m_sim_ptr->repo();
-    const auto nlevels = repo.num_active_levels();
-    const auto geom = m_sim_ptr->mesh().Geom();
-
-    // Get blanking for cells
-    const auto& iblank_cell = repo.get_int_field("iblank_cell");
-
-    // Get fields that will be modified or used
-    Field* rho0{nullptr};
-    auto& rho = repo.get_field("density");
-    auto& gp = repo.get_field("gp");
-    if (m_mphase->perturb_pressure()) {
-        rho0 = &(m_sim_ptr->repo().get_field("reference_density"));
-    } else {
-        // Point to existing field, won't be used
-        rho0 = &rho;
-    }
-
-    // Replace initial gp with best guess (hydrostatic)
-    for (int lev = 0; lev < nlevels; ++lev) {
-        overset_ops::replace_gradp_hydrostatic(
-            gp(lev), rho(lev), (*rho0)(lev), iblank_cell(lev),
-            m_mphase->gravity()[2], m_mphase->perturb_pressure());
-    }
-    amrex::Gpu::synchronize();
-}
-
 void OversetOps::replace_masked_gradp()
 {
     const auto& repo = m_sim_ptr->repo();

amr-wind/overset/TiogaInterface.H

@@ -131,8 +131,6 @@ private:
 
     std::vector<std::string> m_cell_vars;
     std::vector<std::string> m_node_vars;
-
-    bool m_disable_nodal_proj{false};
 };
 
 } // namespace amr_wind

amr-wind/overset/TiogaInterface.cpp

@@ -3,63 +3,14 @@
 #include "amr-wind/core/FieldRepo.H"
 #include "amr-wind/equation_systems/PDEBase.H"
 #include "amr-wind/core/field_ops.H"
+#include "amr-wind/overset/overset_ops_routines.H"
 #include "amr-wind/utilities/IOManager.H"
 #include "AMReX_ParmParse.H"
 
 #include <memory>
 #include <numeric>
 namespace amr_wind {
 
-namespace {
-
-/** Convert iblanks to AMReX mask
- *
- *  \f{align}
- *  \mathrm{mask}_{i,j,k} = \begin{cases}
- *  1 & \mathrm{IBLANK}_{i, j, k} = 1 \\
- *  0 & \mathrm{IBLANK}_{i, j, k} \leq 0
- *  \end{cases}
- *  \f}
- */
-void iblank_to_mask(const IntField& iblank, IntField& maskf)
-{
-    const auto& nlevels = iblank.repo().mesh().finestLevel() + 1;
-
-    for (int lev = 0; lev < nlevels; ++lev) {
-        const auto& ibl = iblank(lev);
-        auto& mask = maskf(lev);
-
-        const auto& ibarrs = ibl.const_arrays();
-        const auto& marrs = mask.arrays();
-        amrex::ParallelFor(
-            ibl, ibl.n_grow,
-            [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
-                marrs[nbx](i, j, k) = amrex::max(ibarrs[nbx](i, j, k), 0);
-            });
-    }
-    amrex::Gpu::synchronize();
-}
-
-void iblank_to_mask_hole(const IntField& iblank, IntField& maskf)
-{
-    const auto& nlevels = iblank.repo().mesh().finestLevel() + 1;
-
-    for (int lev = 0; lev < nlevels; ++lev) {
-        const auto& ibl = iblank(lev);
-        auto& mask = maskf(lev);
-
-        const auto& ibarrs = ibl.const_arrays();
-        const auto& marrs = mask.arrays();
-        amrex::ParallelFor(
-            ibl, ibl.n_grow,
-            [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
-                marrs[nbx](i, j, k) = std::abs(ibarrs[nbx](i, j, k));
-            });
-    }
-    amrex::Gpu::synchronize();
-}
-} // namespace
-
 AMROversetInfo::AMROversetInfo(const int nglobal, const int nlocal)
     : level(nglobal)
     , mpi_rank(nglobal)
@@ -93,9 +44,6 @@ TiogaInterface::TiogaInterface(CFDSim& sim)
           FieldLoc::NODE))
 {
     m_sim.io_manager().register_output_int_var(m_iblank_cell.name());
-
-    amrex::ParmParse pp("Overset");
-    pp.query("disable_coupled_nodal_proj", m_disable_nodal_proj);
 }
 
 // clang-format on
@@ -151,12 +99,8 @@ void TiogaInterface::post_overset_conn_work()
         m_iblank_node(lev).FillBoundary(m_sim.mesh().Geom()[lev].periodicity());
     }
 
-    iblank_to_mask(m_iblank_cell, m_mask_cell);
-    if (m_disable_nodal_proj) {
-        iblank_to_mask_hole(m_iblank_node, m_mask_node);
-    } else {
-        iblank_to_mask(m_iblank_node, m_mask_node);
-    }
+    overset_ops::iblank_to_mask(m_iblank_cell, m_mask_cell);
+    overset_ops::iblank_to_mask(m_iblank_node, m_mask_node);
 
     // Update equation systems after a connectivity update
     m_sim.pde_manager().icns().post_regrid_actions();