Multiphase Hybrid - FY24 Q2 milestone (sharpening the pressure field)

amr-wind/overset/CMakeLists.txt

@@ -2,4 +2,5 @@ target_sources(${amr_wind_lib_name}
   PRIVATE
   TiogaInterface.cpp
   OversetOps.cpp
+  overset_ops_routines.cpp
   )

amr-wind/overset/OversetOps.H

@@ -21,17 +21,15 @@ private:
     // Functions called within public functions
     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};
-    // 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};
+    // To avoid using sharpened pressure, use hydrostatic pressure
+    bool m_use_hydrostatic_gradp{false};
 
     // Coupling options
     bool m_disable_nodal_proj{false};

amr-wind/overset/OversetOps.cpp

@@ -22,19 +22,33 @@ 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);
 
-    amrex::ParmParse pp_icns("ICNS");
-    pp_icns.query("use_perturb_pressure", m_perturb_p);
-
     m_vof_exists = (*m_sim_ptr).repo().field_exists("vof");
     if (m_vof_exists) {
         m_mphase = &(*m_sim_ptr).physics_manager().get<MultiPhase>();
+
+        // Check combination of pressure options
+        if (m_mphase->perturb_pressure() &&
+            !m_mphase->reconstruct_true_pressure()) {
+            amrex::Abort(
+                "OversetOps: perturbational pressure is turned on, but true "
+                "pressure reconstruction is turned off. This approach will be "
+                "incorrect when coupling with Nalu-Wind.");
+        }
+        if (!m_mphase->perturb_pressure() &&
+            m_mphase->reconstruct_true_pressure()) {
+            amrex::Print()
+                << "WARNING (OversetOps): true pressure reconstruction is "
+                   "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);
@@ -45,7 +59,6 @@ void OversetOps::initialize(CFDSim& sim)
 
 void OversetOps::pre_advance_work()
 {
-    // Pressure gradient not updated for current multiphase approach
     if (!(m_vof_exists && m_use_hydrostatic_gradp)) {
         // Update pressure gradient using updated overset pressure field
         update_gradp();
@@ -57,12 +70,19 @@ void OversetOps::pre_advance_work()
         if (m_use_hydrostatic_gradp) {
             // Use hydrostatic pressure gradient
             set_hydrostatic_gradp();
+        } else {
+            if (m_mphase->perturb_pressure()) {
+                // Modify to be consistent with internal source terms
+                form_perturb_pressure();
+            }
+            // Update pressure gradient using sharpened pressure field
+            update_gradp();
         }
     }
 
     // If pressure gradient will be replaced, store current pressure gradient
     if (m_replace_gradp_postsolve) {
-        auto& gp = (*m_sim_ptr).repo().get_field("gp");
+        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(
@@ -154,6 +174,11 @@ void OversetOps::parameter_output() const
                        << "---- 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
@@ -190,29 +215,50 @@ void OversetOps::sharpen_nalu_data()
     auto& levelset = repo.get_field("levelset");
     auto& rho = repo.get_field("density");
     auto& velocity = repo.get_field("velocity");
+    auto& gp_noghost = repo.get_field("gp");
+    auto& p = repo.get_field("p");
 
-    // Create scratch fields for fluxes
-    // 5 components are vof, density, and 3 of velocity
-    auto flux_x = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::XFACE);
-    auto flux_y = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::YFACE);
-    auto flux_z = repo.create_scratch_field(5, 0, amr_wind::FieldLoc::ZFACE);
-    // Create scratch field for approximate signed distance function and grad
-    // (components 0-2 are gradient, 3 is asdf)
-    auto normal_vec = repo.create_scratch_field(3, vof.num_grow()[0] - 1);
+    // 9 components are vof, density, 3 of velocity, 3 of gp, and psource flag
+    auto flux_x = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::XFACE);
+    auto flux_y = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::YFACE);
+    auto flux_z = repo.create_scratch_field(9, 1, amr_wind::FieldLoc::ZFACE);
 
+    auto p_src = repo.create_scratch_field(1, 0, amr_wind::FieldLoc::NODE);
+    auto normal_vec = repo.create_scratch_field(3, vof.num_grow()[0] - 1);
     auto target_vof = repo.create_scratch_field(1, vof.num_grow()[0]);
 
-    // Create levelset field
+    // Sharpening fluxes (at faces) have 1 ghost, requiring fields to have >= 2
+    auto gp_scr = repo.create_scratch_field(3, 2);
+    auto& gp = *gp_scr;
+
+    // Give initial max possible value of pseudo-velocity scale
+    const auto dx_lev0 = (geom[0]).CellSizeArray();
+    const amrex::Real max_pvscale =
+        std::min(std::min(dx_lev0[0], dx_lev0[1]), dx_lev0[2]);
+    amrex::Real pvscale = max_pvscale;
+
+    // Prep things that do not change with iterations
     for (int lev = 0; lev < nlevels; ++lev) {
         // Thickness used here is user parameter, whatever works best
-        auto dx = (geom[lev]).CellSizeArray();
+        const auto dx = (geom[lev]).CellSizeArray();
         const amrex::Real i_th =
             m_relative_length_scale * std::cbrt(dx[0] * dx[1] * dx[2]);
 
         // Populate approximate signed distance function
         overset_ops::populate_psi(levelset(lev), vof(lev), i_th, m_asdf_tiny);
+
+        gp(lev).setVal(0.0);
+        amrex::MultiFab::Copy(gp(lev), gp_noghost(lev), 0, 0, 3, 0);
+        gp(lev).FillBoundary(geom[lev].periodicity());
+
+        // Get pseudo-velocity scale, proportional to smallest dx in iblank
+        const amrex::Real pvscale_lev =
+            overset_ops::calculate_pseudo_velocity_scale(
+                iblank_cell(lev), dx, max_pvscale);
+        pvscale = std::min(pvscale, pvscale_lev);
     }
     amrex::Gpu::synchronize();
+    amrex::ParallelDescriptor::ReduceRealMin(pvscale);
 
     // Convert levelset to vof to get target_vof
     m_mphase->levelset2vof(iblank_cell, *target_vof);
@@ -262,13 +308,14 @@ void OversetOps::sharpen_nalu_data()
             // Sharpening fluxes for vof, density, and momentum
             overset_ops::populate_sharpen_fluxes(
                 (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
-                (*target_vof)(lev), (*normal_vec)(lev), velocity(lev),
-                m_upw_margin, m_mphase->rho1(), m_mphase->rho2());
+                (*target_vof)(lev), (*normal_vec)(lev), velocity(lev), gp(lev),
+                rho(lev), pvscale, m_upw_margin, m_mphase->rho1(),
+                m_mphase->rho2());
 
             // Process fluxes
-            overset_ops::process_fluxes(
-                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev),
-                iblank_cell(lev));
+            overset_ops::process_fluxes_calc_src(
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev),
+                vof(lev), iblank_cell(lev));
 
             // Measure convergence to determine if loop can stop
             if (calc_convg) {
@@ -282,11 +329,9 @@ void OversetOps::sharpen_nalu_data()
 
         // Average down fluxes across levels for consistency
         for (int lev = nlevels - 1; lev > 0; --lev) {
-            amrex::IntVect rr =
-                geom[lev].Domain().size() / geom[lev - 1].Domain().size();
             amrex::average_down_faces(
-                GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1], rr,
-                geom[lev - 1]);
+                GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1],
+                repo.mesh().refRatio(lev), geom[lev - 1]);
         }
 
         // Get pseudo dt (dtau)
@@ -305,15 +350,19 @@ void OversetOps::sharpen_nalu_data()
 
         // Apply fluxes
         for (int lev = 0; lev < nlevels; ++lev) {
+            const auto dx = (geom[lev]).CellSizeArray();
+
             overset_ops::apply_fluxes(
-                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), vof(lev),
-                rho(lev), velocity(lev), ptfac, m_vof_tol);
+                (*flux_x)(lev), (*flux_y)(lev), (*flux_z)(lev), (*p_src)(lev),
+                vof(lev), rho(lev), velocity(lev), gp(lev), p(lev), dx, ptfac,
+                m_vof_tol);
+
+            vof(lev).FillBoundary(geom[lev].periodicity());
+            velocity(lev).FillBoundary(geom[lev].periodicity());
+            gp(lev).FillBoundary(geom[lev].periodicity());
         }
-        amrex::Gpu::synchronize();
 
-        // Fillpatch for ghost cells
-        vof.fillpatch((*m_sim_ptr).time().current_time());
-        velocity.fillpatch((*m_sim_ptr).time().current_time());
+        amrex::Gpu::synchronize();
 
         // Update density (fillpatch built in)
         m_mphase->set_density_via_vof();
@@ -329,6 +378,9 @@ void OversetOps::sharpen_nalu_data()
         }
     }
 
+    // Fillpatch for pressure to make sure pressure stencil has all points
+    p.fillpatch((*m_sim_ptr).time().current_time());
+
     // Purely for debugging via visualization, should be removed later
     // Currently set up to overwrite the levelset field (not used as time
     // evolves) with the post-sharpening velocity magnitude
@@ -338,6 +390,16 @@ void OversetOps::sharpen_nalu_data()
     amrex::Gpu::synchronize();
 }
 
+void OversetOps::form_perturb_pressure()
+{
+    auto& pressure = (*m_sim_ptr).repo().get_field("p");
+    const auto& p0 = (*m_sim_ptr).repo().get_field("reference_pressure");
+    for (int lev = 0; lev < (*m_sim_ptr).repo().num_active_levels(); lev++) {
+        amrex::MultiFab::Subtract(
+            pressure(lev), p0(lev), 0, 0, 1, pressure.num_grow()[0]);
+    }
+}
+
 void OversetOps::set_hydrostatic_gradp()
 {
     const auto& repo = (*m_sim_ptr).repo();
@@ -351,7 +413,7 @@ void OversetOps::set_hydrostatic_gradp()
     Field* rho0{nullptr};
     auto& rho = repo.get_field("density");
     auto& gp = repo.get_field("gp");
-    if (m_perturb_p) {
+    if (m_mphase->perturb_pressure()) {
         rho0 = &((*m_sim_ptr).repo().get_field("reference_density"));
     } else {
         // Point to existing field, won't be used
@@ -362,7 +424,7 @@ void OversetOps::set_hydrostatic_gradp()
     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_perturb_p);
+            m_mphase->gravity()[2], m_mphase->perturb_pressure());
     }
 }