More multiphase overset improvements (#1182)

* proper use of refRatio
* remove pressure initialization from sloshing tank
* correct approach on first timestep
* removing unused function declarations

Author: mbkuhn

amr-wind/incflo.H

@@ -143,11 +143,6 @@ public:
         amrex::Real scaling_factor,
         bool incremental);
 
-    void UpdateGradP(
-        amrex::Vector<amrex::MultiFab const*> density,
-        amrex::Real time,
-        amrex::Real scaling_factor);
-
     //! Initialize Physics instances as well as PDEs (include turbulence models)
     void init_physics_and_pde();
 
@@ -204,7 +199,6 @@ private:
     //
     ///////////////////////////////////////////////////////////////////////////
 
-    void set_background_pressure();
     void ReadParameters();
     void InitialProjection();
     void InitialIterations();

amr-wind/overset/OversetOps.cpp

@@ -30,9 +30,9 @@ void OversetOps::initialize(CFDSim& sim)
 
     pp.query("verbose", m_verbose);
 
-    m_vof_exists = (*m_sim_ptr).repo().field_exists("vof");
+    m_vof_exists = m_sim_ptr->repo().field_exists("vof");
     if (m_vof_exists) {
-        m_mphase = &(*m_sim_ptr).physics_manager().get<MultiPhase>();
+        m_mphase = &m_sim_ptr->physics_manager().get<MultiPhase>();
 
         // Check combination of pressure options
         if (m_mphase->perturb_pressure() &&
@@ -51,7 +51,7 @@ void OversetOps::initialize(CFDSim& sim)
         }
     }
     if (m_replace_gradp_postsolve) {
-        m_gp_copy = &(*m_sim_ptr).repo().declare_field("gp_copy", 3);
+        m_gp_copy = &m_sim_ptr->repo().declare_field("gp_copy", 3);
     }
 
     parameter_output();
@@ -60,6 +60,14 @@ void OversetOps::initialize(CFDSim& sim)
 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();
+            }
+        }
         // Update pressure gradient using updated overset pressure field
         update_gradp();
     }
@@ -71,20 +79,15 @@ void OversetOps::pre_advance_work()
             // 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) {
-        const auto& gp = (*m_sim_ptr).repo().get_field("gp");
-        for (int lev = 0; lev < (*m_sim_ptr).repo().num_active_levels();
-             ++lev) {
+        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());
@@ -97,30 +100,30 @@ void OversetOps::update_gradp()
     BL_PROFILE("amr-wind::OversetOps::update_gradp");
 
     // Get relevant fields
-    auto& grad_p = (*m_sim_ptr).repo().get_field("gp");
-    auto& pressure = (*m_sim_ptr).repo().get_field("p");
-    auto& velocity = (*m_sim_ptr).repo().get_field("velocity");
+    auto& grad_p = m_sim_ptr->repo().get_field("gp");
+    auto& pressure = m_sim_ptr->repo().get_field("p");
+    auto& velocity = m_sim_ptr->repo().get_field("velocity");
 
     // Set up projection object
     std::unique_ptr<Hydro::NodalProjector> nodal_projector;
     // Boundary conditions from field, periodicity
     auto bclo = nodal_projection::get_projection_bc(
         amrex::Orientation::low, pressure,
-        (*m_sim_ptr).mesh().Geom(0).isPeriodic());
+        m_sim_ptr->mesh().Geom(0).isPeriodic());
     auto bchi = nodal_projection::get_projection_bc(
         amrex::Orientation::high, pressure,
-        (*m_sim_ptr).mesh().Geom(0).isPeriodic());
+        m_sim_ptr->mesh().Geom(0).isPeriodic());
     // Velocity multifab is needed for proper initialization, but only the size
     // matters for the purpose of calculating gradp, the values do not matter
-    int finest_level = (*m_sim_ptr).mesh().finestLevel();
+    int finest_level = m_sim_ptr->mesh().finestLevel();
     Vector<MultiFab*> vel;
     for (int lev = 0; lev <= finest_level; ++lev) {
         vel.push_back(&(velocity(lev)));
     }
     amr_wind::MLMGOptions options("nodal_proj");
     // Create nodal projector with unity scaling factor for simplicity
     nodal_projector = std::make_unique<Hydro::NodalProjector>(
-        vel, 1.0, (*m_sim_ptr).mesh().Geom(0, finest_level), options.lpinfo());
+        vel, 1.0, m_sim_ptr->mesh().Geom(0, finest_level), options.lpinfo());
     // Set MLMG and NodalProjector options
     options(*nodal_projector);
     nodal_projector->setDomainBC(bclo, bchi);
@@ -203,9 +206,9 @@ void OversetOps::parameter_output() const
 
 void OversetOps::sharpen_nalu_data()
 {
-    const auto& repo = (*m_sim_ptr).repo();
+    const auto& repo = m_sim_ptr->repo();
     const auto nlevels = repo.num_active_levels();
-    const auto geom = (*m_sim_ptr).mesh().Geom();
+    const auto geom = m_sim_ptr->mesh().Geom();
 
     // Get blanking for cells
     const auto& iblank_cell = repo.get_int_field("iblank_cell");
@@ -333,7 +336,7 @@ void OversetOps::sharpen_nalu_data()
         for (int lev = nlevels - 1; lev > 0; --lev) {
             amrex::average_down_faces(
                 GetArrOfConstPtrs(fluxes[lev]), fluxes[lev - 1],
-                repo.mesh().refRatio(lev), geom[lev - 1]);
+                repo.mesh().refRatio(lev - 1), geom[lev - 1]);
         }
 
         // Get pseudo dt (dtau)
@@ -382,7 +385,7 @@ void OversetOps::sharpen_nalu_data()
     }
 
     // Fillpatch for pressure to make sure pressure stencil has all points
-    p.fillpatch((*m_sim_ptr).time().current_time());
+    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
@@ -395,19 +398,19 @@ void OversetOps::sharpen_nalu_data()
 
 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++) {
+    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();
+    const auto& repo = m_sim_ptr->repo();
     const auto nlevels = repo.num_active_levels();
-    const auto geom = (*m_sim_ptr).mesh().Geom();
+    const auto geom = m_sim_ptr->mesh().Geom();
 
     // Get blanking for cells
     const auto& iblank_cell = repo.get_int_field("iblank_cell");
@@ -417,7 +420,7 @@ void OversetOps::set_hydrostatic_gradp()
     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"));
+        rho0 = &(m_sim_ptr->repo().get_field("reference_density"));
     } else {
         // Point to existing field, won't be used
         rho0 = &rho;
@@ -433,11 +436,11 @@ void OversetOps::set_hydrostatic_gradp()
 
 void OversetOps::replace_masked_gradp()
 {
-    const auto& repo = (*m_sim_ptr).repo();
+    const auto& repo = m_sim_ptr->repo();
     const auto nlevels = repo.num_active_levels();
 
     // Get timestep
-    const amrex::Real dt = (*m_sim_ptr).time().delta_t();
+    const amrex::Real dt = m_sim_ptr->time().delta_t();
 
     // Get blanking for cells
     const auto& iblank_cell = repo.get_int_field("iblank_cell");

amr-wind/physics/multiphase/SloshingTank.H

@@ -35,7 +35,6 @@ public:
 private:
     Field& m_velocity;
     Field& m_levelset;
-    Field& m_pressure;
 
     //! Initial free surface amplitude magnitude
     amrex::Real m_amplitude{0.1};
@@ -45,14 +44,6 @@ private:
 
     //! Initial zero-level free-surface water depth
     amrex::Real m_waterlevel{0.0};
-
-    //! Stuff to get from MultiPhase physics
-    amrex::Real m_rho1{std::numeric_limits<amrex::Real>::quiet_NaN()};
-    amrex::Real m_rho2{std::numeric_limits<amrex::Real>::quiet_NaN()};
-    amrex::Vector<amrex::Real> m_gravity{
-        std::numeric_limits<amrex::Real>::quiet_NaN(),
-        std::numeric_limits<amrex::Real>::quiet_NaN(),
-        std::numeric_limits<amrex::Real>::quiet_NaN()};
 };
 
 } // namespace amr_wind

amr-wind/physics/multiphase/SloshingTank.cpp

@@ -1,5 +1,4 @@
 #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"
@@ -10,30 +9,22 @@ 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);
-
-    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
- *  simulation. Initializing pressure has little effect for pure AMR-Wind
- *  simulations, but it improves initialization in the overset solver.
+ *  simulation.
  */
 void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
 {
     auto& velocity = m_velocity(level);
     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();
@@ -42,14 +33,10 @@ 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];
 
     for (amrex::MFIter mfi(levelset); mfi.isValid(); ++mfi) {
         const auto& vbx = mfi.validbox();
         auto phi = levelset.array(mfi);
-        auto p = pressure.array(mfi);
 
         amrex::ParallelFor(
             vbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
@@ -63,29 +50,6 @@ void SloshingTank::initialize_fields(int level, const amrex::Geometry& geom)
                                         std::pow(y - problo[1] - 0.5 * Ly, 2)));
                 phi(i, j, k) = z0 - z;
             });
-        amrex::Box const& nbx = mfi.grownnodaltilebox();
-        amrex::ParallelFor(
-            nbx, [=] AMREX_GPU_DEVICE(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
-                const amrex::Real ih_g =
-                    amrex::max(0.0, amrex::min(probhi[2] - z0, probhi[2] - z));
-                const 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(i, j, k) = -irho * grav_z;
-            });
     }
 }