Initialize pressure for sloshing tank, waves cases (overset purposes) (#1687)

* make pressure initialization available

* tweaks to interface_band

* incorporate iblank_to_mask_vof

* zero velocities inside objects

* make vof-dependent node mask when nodal_proj coupling not disabled

* correcting arguments

* getting pressure changes back in

* behavior matches now with perturb pressure

* formatting

* remove some unnecessary changes

* reinsert sync

* turn off sloshing tank pressure by default

Author: mbkuhn

amr-wind/overset/OversetOps.cpp

@@ -56,20 +56,16 @@ void OversetOps::initialize(CFDSim& sim)
 
 void OversetOps::pre_advance_work()
 {
-    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();
 
     if (m_vof_exists) {
         // Reinitialize fields
         sharpen_nalu_data();
+        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();
         // Calculate vof-dependent node mask

amr-wind/physics/multiphase/MultiPhase.H

@@ -104,6 +104,8 @@ private:
     // Reconstructing true pressure field at end of timestep
     bool m_reconstruct_true_pressure{false};
 
+    bool m_init_p{false};
+
     // Info to create rho0
     amrex::Real m_water_level0{0.0};
     // Info to reconstruct true pressure

amr-wind/physics/multiphase/MultiPhase.cpp

@@ -73,6 +73,9 @@ MultiPhase::MultiPhase(CFDSim& sim)
         amrex::Print() << "WARNING: single-phase density has been specified "
                           "but will not be used! (MultiPhase physics)\n";
     }
+
+    pp_multiphase.query("initialize_pressure", m_init_p);
+
     // Always populate gravity
     pp_incflo.queryarr("gravity", m_gravity);
 }
@@ -131,6 +134,20 @@ void MultiPhase::post_init_actions()
                 m_sim.mesh().Geom());
         }
     }
+
+    if (m_init_p && !is_wlev) {
+        amrex::Abort(
+            "Initialize pressure requested, but physics case does not "
+            "specify water level.");
+    }
+    // Make p field if both are specified
+    if (m_init_p && is_wlev) {
+        // Initialize rho0 field for perturbational density, pressure
+        auto& p = m_sim.repo().get_field("p");
+        hydrostatic::define_p0(
+            p, m_rho1, m_rho2, m_water_level0, m_gravity[2],
+            m_sim.mesh().Geom());
+    }
 }
 
 void MultiPhase::post_regrid_actions()

amr-wind/physics/multiphase/SloshingTank.H

@@ -35,6 +35,7 @@ public:
 private:
     Field& m_velocity;
     Field& m_levelset;
+    Field& m_pressure;
 
     //! Initial free surface amplitude magnitude
     amrex::Real m_amplitude{0.1};
@@ -44,6 +45,13 @@ private:
 
     //! Initial zero-level free-surface water depth
     amrex::Real m_waterlevel{0.0};
+
+    bool m_init_p{false};
+
+    //! Stuff to get from MultiPhase physics
+    amrex::Real m_rho1{1000.};
+    amrex::Real m_rho2{1.};
+    amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
 };
 
 } // namespace amr_wind

amr-wind/physics/multiphase/SloshingTank.cpp

@@ -1,4 +1,5 @@
 #include "amr-wind/physics/multiphase/SloshingTank.H"
+#include "amr-wind/physics/multiphase/MultiPhase.H"
 #include "amr-wind/CFDSim.H"
 #include "AMReX_ParmParse.H"
 #include "amr-wind/fvm/gradient.H"
@@ -9,11 +10,19 @@ namespace amr_wind {
 SloshingTank::SloshingTank(CFDSim& sim)
     : m_velocity(sim.repo().get_field("velocity"))
     , m_levelset(sim.repo().get_field("levelset"))
+    , m_pressure(sim.repo().get_field("p"))
 {
     amrex::ParmParse pp(identifier());
     pp.query("amplitude", m_amplitude);
     pp.query("peak_enhance", m_kappa);
     pp.query("water_level", m_waterlevel);
+
+    pp.query("initialize_pressure", m_init_p);
+
+    const auto& mphase = sim.physics_manager().get<MultiPhase>();
+    m_rho1 = mphase.rho1();
+    m_rho2 = mphase.rho2();
+    m_gravity = mphase.gravity();
 }
 
 /** Initialize the velocity and levelset fields at the beginning of the
@@ -25,6 +34,7 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
     velocity.setVal(0.0);
 
     auto& levelset = m_levelset(level);
+    auto& pressure = m_pressure(level);
     const auto& dx = geom.CellSizeArray();
     const auto& problo = geom.ProbLoArray();
     const auto& probhi = geom.ProbHiArray();
@@ -33,7 +43,12 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
     const amrex::Real water_level = m_waterlevel;
     const amrex::Real Lx = probhi[0] - problo[0];
     const amrex::Real Ly = probhi[1] - problo[1];
+    const amrex::Real rho1 = m_rho1;
+    const amrex::Real rho2 = m_rho2;
+    const amrex::Real grav_z = m_gravity[2];
+
     const auto& phi_arrs = levelset.arrays();
+    const auto& p = pressure.arrays();
 
     amrex::ParallelFor(
         levelset, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
@@ -47,6 +62,32 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
                                     std::pow(y - problo[1] - 0.5 * Ly, 2)));
             phi_arrs[nbx](i, j, k) = z0 - z;
         });
+
+    if (m_init_p) {
+        amrex::ParallelFor(
+            pressure, amrex::IntVect(0),
+            [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
+                // For pressure nodes, no offset
+                const amrex::Real x = problo[0] + i * dx[0];
+                const amrex::Real y = problo[1] + j * dx[1];
+                const amrex::Real z = problo[2] + k * dx[2];
+                const amrex::Real z0 =
+                    water_level +
+                    Amp * std::exp(
+                              -kappa * (std::pow(x - problo[0] - 0.5 * Lx, 2) +
+                                        std::pow(y - problo[1] - 0.5 * Ly, 2)));
+                // Integrated (top-down in z) phase heights to pressure node
+                amrex::Real ih_g =
+                    amrex::max(0.0, amrex::min(probhi[2] - z0, probhi[2] - z));
+                amrex::Real ih_l =
+                    amrex::max(0.0, amrex::min(z0 - z, z0 - problo[2]));
+                // Integrated rho at pressure node
+                const amrex::Real irho = rho1 * ih_l + rho2 * ih_g;
+
+                // Add term to reference pressure
+                p[nbx](i, j, k) = -irho * grav_z;
+            });
+    }
     amrex::Gpu::streamSynchronize();
 }