From 6ef259a3eb60cf1bc80867e544b6cf56bc228fd7 Mon Sep 17 00:00:00 2001 From: Michael B Kuhn <31661049+mbkuhn@users.noreply.github.com> Date: Wed, 30 Aug 2023 14:04:39 -0600 Subject: [PATCH] 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 | 5 + amr-wind/incflo_advance.cpp | 7 ++ .../incflo_apply_nodal_projection.cpp | 116 ++++++++++++++++++ submods/AMReX-Hydro | 2 +- 4 files changed, 129 insertions(+), 1 deletion(-) diff --git a/amr-wind/incflo.H b/amr-wind/incflo.H index 3105b2f37d..3eca0d7dff 100644 --- a/amr-wind/incflo.H +++ b/amr-wind/incflo.H @@ -145,6 +145,11 @@ public: amrex::Real scaling_factor, bool incremental); + void UpdateGradP( + amrex::Vector density, + amrex::Real time, + amrex::Real scaling_factor); + //! Initialize Physics instances as well as PDEs (include turbulence models) void init_physics_and_pde(); diff --git a/amr-wind/incflo_advance.cpp b/amr-wind/incflo_advance.cpp index a230c7a7e0..b4ef5d95ac 100644 --- a/amr-wind/incflo_advance.cpp +++ b/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 // ************************************************************************************* diff --git a/amr-wind/projection/incflo_apply_nodal_projection.cpp b/amr-wind/projection/incflo_apply_nodal_projection.cpp index b7c1a9c5a2..c042b7e638 100644 --- a/amr-wind/projection/incflo_apply_nodal_projection.cpp +++ b/amr-wind/projection/incflo_apply_nodal_projection.cpp @@ -435,3 +435,119 @@ void incflo::ApplyProjection( } } } + +void incflo::UpdateGradP( + Vector 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 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 const& sig = sigma[lev].array(mfi); + Array4 const& rho = density[lev]->const_array(mfi); + amrex::Array4 fac = + mesh_mapping ? ((*mesh_fac)(lev).const_array(mfi)) + : amrex::Array4(); + amrex::Array4 detJ = + mesh_mapping ? ((*mesh_detJ)(lev).const_array(mfi)) + : amrex::Array4(); + + 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 nodal_projector; + + auto bclo = get_projection_bc(Orientation::low); + auto bchi = get_projection_bc(Orientation::high); + + Vector 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( + 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( + 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 const& gp_lev = grad_p(lev).array(mfi); + Array4 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 +} diff --git a/submods/AMReX-Hydro b/submods/AMReX-Hydro index bc4e5f64c9..c6fed55091 160000 --- a/submods/AMReX-Hydro +++ b/submods/AMReX-Hydro @@ -1 +1 @@ -Subproject commit bc4e5f64c9ac3072793ebf426b422be91231652c +Subproject commit c6fed550913cc477cd964e8d137502da5555090c