overset: calc pressure gradient before proceeding (#901)

- update AMReX-Hydro for new function
- use new function to get nalu pressure into gp

amr-wind/incflo.H

@@ -145,6 +145,11 @@ 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();
 

amr-wind/incflo_advance.cpp

@@ -190,6 +190,13 @@ 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
+    if (sim().has_overset()) {
+        UpdateGradP(
+            (density_old).vec_const_ptrs(), m_time.current_time(),
+            m_time.deltaT());
+    }
+
     // *************************************************************************************
     // Compute viscosity / diffusive coefficients
     // *************************************************************************************

amr-wind/projection/incflo_apply_nodal_projection.cpp

@@ -435,3 +435,119 @@ void incflo::ApplyProjection(
         }
     }
 }
+
+void incflo::UpdateGradP(
+    Vector<MultiFab const*> density, Real time, Real scaling_factor)
+{
+    BL_PROFILE("amr-wind::incflo::UpdateGradP");
+
+    bool variable_density =
+        (!m_sim.pde_manager().constant_density() ||
+         m_sim.physics_manager().contains("MultiPhase"));
+
+    bool mesh_mapping = m_sim.has_mesh_mapping();
+
+    auto& grad_p = m_repo.get_field("gp");
+    auto& pressure = m_repo.get_field("p");
+    auto& velocity = icns().fields().field;
+    amr_wind::Field const* mesh_fac =
+        mesh_mapping
+            ? &(m_repo.get_mesh_mapping_field(amr_wind::FieldLoc::CELL))
+            : nullptr;
+    amr_wind::Field const* mesh_detJ =
+        mesh_mapping ? &(m_repo.get_mesh_mapping_detJ(amr_wind::FieldLoc::CELL))
+                     : nullptr;
+
+    // Create sigma while accounting for mesh mapping
+    // sigma = 1/(fac^2)*J * dt/rho
+    Vector<amrex::MultiFab> sigma(finest_level + 1);
+    if (variable_density || mesh_mapping) {
+        int ncomp = mesh_mapping ? AMREX_SPACEDIM : 1;
+        for (int lev = 0; lev <= finest_level; ++lev) {
+            sigma[lev].define(
+                grids[lev], dmap[lev], ncomp, 0, MFInfo(), Factory(lev));
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+            for (MFIter mfi(sigma[lev], TilingIfNotGPU()); mfi.isValid();
+                 ++mfi) {
+                Box const& bx = mfi.tilebox();
+                Array4<Real> const& sig = sigma[lev].array(mfi);
+                Array4<Real const> const& rho = density[lev]->const_array(mfi);
+                amrex::Array4<amrex::Real const> fac =
+                    mesh_mapping ? ((*mesh_fac)(lev).const_array(mfi))
+                                 : amrex::Array4<amrex::Real const>();
+                amrex::Array4<amrex::Real const> detJ =
+                    mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi))
+                                 : amrex::Array4<amrex::Real const>();
+
+                amrex::ParallelFor(
+                    bx, ncomp,
+                    [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
+                        amrex::Real fac_cc =
+                            mesh_mapping ? (fac(i, j, k, n)) : 1.0;
+                        amrex::Real det_j =
+                            mesh_mapping ? (detJ(i, j, k)) : 1.0;
+                        sig(i, j, k, n) = std::pow(fac_cc, -2.) * det_j *
+                                          scaling_factor / rho(i, j, k);
+                    });
+            }
+        }
+    }
+
+    // Perform projection
+    std::unique_ptr<Hydro::NodalProjector> nodal_projector;
+
+    auto bclo = get_projection_bc(Orientation::low);
+    auto bchi = get_projection_bc(Orientation::high);
+
+    Vector<MultiFab*> vel;
+    for (int lev = 0; lev <= finest_level; ++lev) {
+        vel.push_back(&(velocity(lev)));
+        vel[lev]->setBndry(0.0);
+        set_inflow_velocity(lev, time, *vel[lev], 1);
+    }
+
+    amr_wind::MLMGOptions options("nodal_proj");
+
+    if (variable_density || mesh_mapping) {
+        nodal_projector = std::make_unique<Hydro::NodalProjector>(
+            vel, GetVecOfConstPtrs(sigma), Geom(0, finest_level),
+            options.lpinfo());
+    } else {
+        amrex::Real rho_0 = 1.0;
+        amrex::ParmParse pp("incflo");
+        pp.query("density", rho_0);
+
+        nodal_projector = std::make_unique<Hydro::NodalProjector>(
+            vel, scaling_factor / rho_0, Geom(0, finest_level),
+            options.lpinfo());
+    }
+
+    // Set MLMG and NodalProjector options
+    options(*nodal_projector);
+    nodal_projector->setDomainBC(bclo, bchi);
+
+    // Recalculate gradphi with fluxes
+    auto gradphi = nodal_projector->calcGradPhi(pressure.vec_ptrs());
+
+    // Transfer pressure gradient to gp field
+    for (int lev = 0; lev <= finest_level; lev++) {
+
+#ifdef AMREX_USE_OMP
+#pragma omp parallel if (Gpu::notInLaunchRegion())
+#endif
+        for (MFIter mfi(grad_p(lev), TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+            Box const& tbx = mfi.tilebox();
+            Array4<Real> const& gp_lev = grad_p(lev).array(mfi);
+            Array4<Real const> const& gp_proj = gradphi[lev]->const_array(mfi);
+            amrex::ParallelFor(
+                tbx, AMREX_SPACEDIM,
+                [=] AMREX_GPU_DEVICE(int i, int j, int k, int n) noexcept {
+                    gp_lev(i, j, k, n) = gp_proj(i, j, k, n);
+                });
+        }
+    }
+
+    // Average down is unnecessary because it is build into calcGradPhi
+}