Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use multiply_rho true for potential temperature evolution #1498

Draft
wants to merge 8 commits into
base: main
Choose a base branch
from
Draft

Use multiply_rho true for potential temperature evolution #1498

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
e6eaab5
Use multiply_rho true for potential temperature evolution
marchdf Feb 20, 2025
9cc955e
use hybrid_fluxes for scalars, thereby for temperature
mbkuhn Feb 21, 2025
80c60bc
meaningless tweaks
mbkuhn Feb 21, 2025
88e3597
put scalar advection in the loop so that vof post-solve happens befor…
mbkuhn Feb 21, 2025
7ef44fc
remove accidental print statement
mbkuhn Feb 21, 2025
c4441dd
improve a conditional and get source term treatment correct
mbkuhn Feb 24, 2025
c857560
async arena
marchdf Feb 24, 2025
76e2ae7
fix warn
marchdf Feb 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 54 additions & 8 deletions amr-wind/equation_systems/AdvOp_Godunov.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@
#include "AMReX_MultiFabUtil.H"
#include "hydro_utils.H"

#include "amr-wind/equation_systems/vof/vof_momentum_flux.H"

namespace amr_wind::pde {

/** Godunov advection operator for scalar transport equations
Expand Down Expand Up @@ -71,6 +73,16 @@ struct AdvectionOp<
"godunov_type is specified, the default weno_z is used.");
}

// Flux calculation used in multiphase portions of domain
pp.query("mflux_type", mflux_type);
if (amrex::toLower(mflux_type) == "minmod") {
mflux_scheme = godunov::scheme::MINMOD;
} else if (amrex::toLower(mflux_type) == "upwind") {
mflux_scheme = godunov::scheme::UPWIND;
} else {
amrex::Abort("Invalid argument entered for mflux_type.");
}

amrex::ParmParse pp_eq{fields_in.field.base_name()};
pp_eq.query(
"allow_inflow_at_pressure_outflow", m_allow_inflow_on_outflow);
Expand Down Expand Up @@ -115,6 +127,8 @@ struct AdvectionOp<
const auto& den = density.state(fstate);
const auto& den_nph = density.state(amr_wind::FieldState::NPH);

const bool mphase_vof = repo.field_exists("vof");

for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
amrex::MFItInfo mfi_info;
if (amrex::Gpu::notInLaunchRegion()) {
Expand All @@ -135,6 +149,9 @@ struct AdvectionOp<
amrex::Array4<amrex::Real> rhotrac;
amrex::FArrayBox rhotrac_nph_fab;
amrex::Array4<amrex::Real> rhotrac_nph;
amrex::FArrayBox srctrac_over_rho_fab;
amrex::Array4<amrex::Real> srctrac_over_rho;
const auto& src_arr = src_term(lev).const_array(mfi);

if (PDE::multiply_rho) {
auto rhotrac_box =
Expand All @@ -145,6 +162,11 @@ struct AdvectionOp<
rhotrac_nph_fab.resize(
rhotrac_box, PDE::ndim, amrex::The_Async_Arena());
rhotrac_nph = rhotrac_nph_fab.array();
if (mphase_vof) {
srctrac_over_rho_fab.resize(
rhotrac_box, PDE::ndim, amrex::The_Async_Arena());
srctrac_over_rho = srctrac_over_rho_fab.array();
}

amrex::ParallelFor(
rhotrac_box, PDE::ndim,
Expand All @@ -154,6 +176,13 @@ struct AdvectionOp<
rho_arr(i, j, k) * tra_arr(i, j, k, n);
rhotrac_nph(i, j, k, n) =
rho_nph_arr(i, j, k) * tra_nph_arr(i, j, k, n);
if (mphase_vof) {
// Need to remove the effect of the multiply_rho
// routine -- the setup is different for
// mass-consistent multiphase advection
srctrac_over_rho(i, j, k, n) =
src_arr(i, j, k, n) / rho_nph_arr(i, j, k);
}
});
}

Expand Down Expand Up @@ -190,26 +219,41 @@ struct AdvectionOp<

HydroUtils::ComputeFluxesOnBoxFromState(
bx, PDE::ndim, mfi,
(PDE::multiply_rho ? rhotrac : tra_arr),
(PDE::multiply_rho ? rhotrac_nph : tra_nph_arr),
((PDE::multiply_rho && !mphase_vof) ? rhotrac
: tra_arr),
((PDE::multiply_rho && !mphase_vof) ? rhotrac_nph
: tra_nph_arr),
(*flux_x)(lev).array(mfi), (*flux_y)(lev).array(mfi),
(*flux_z)(lev).array(mfi), (*face_x)(lev).array(mfi),
(*face_y)(lev).array(mfi), (*face_z)(lev).array(mfi),
known_edge_state, u_mac(lev).const_array(mfi),
v_mac(lev).const_array(mfi),
w_mac(lev).const_array(mfi), divu,
src_term(lev).const_array(mfi), geom[lev], dt_extrap,
dof_field.bcrec(), dof_field.bcrec_device().data(),
iconserv.data(), godunov_use_ppm,
godunov_use_forces_in_trans, is_velocity,
fluxes_are_area_weighted, advection_type, limiter_type,
m_allow_inflow_on_outflow);
((PDE::multiply_rho && mphase_vof)
? srctrac_over_rho
: src_term(lev).const_array(mfi)),
geom[lev], dt_extrap, dof_field.bcrec(),
dof_field.bcrec_device().data(), iconserv.data(),
godunov_use_ppm, godunov_use_forces_in_trans,
is_velocity, fluxes_are_area_weighted, advection_type,
limiter_type, m_allow_inflow_on_outflow);
} else {
amrex::Abort("Invalid godunov scheme");
}
}
}

// Multiphase flux operations
if (mphase_vof && PDE::multiply_rho) {
// Loop levels
multiphase::hybrid_fluxes(
repo, PDE::ndim, iconserv, (*flux_x), (*flux_y), (*flux_z),
dof_field, dof_nph, src_term, den, den_nph, u_mac, v_mac, w_mac,
dof_field.bcrec(), dof_field.bcrec_device().data(), den.bcrec(),
den.bcrec_device().data(), dt, mflux_scheme,
m_allow_inflow_on_outflow, godunov_use_forces_in_trans, true);
}

amrex::Vector<amrex::Array<amrex::MultiFab*, AMREX_SPACEDIM>> fluxes(
repo.num_active_levels());
for (int lev = 0; lev < repo.num_active_levels(); ++lev) {
Expand Down Expand Up @@ -253,7 +297,9 @@ struct AdvectionOp<
amrex::Gpu::DeviceVector<int> iconserv;

godunov::scheme godunov_scheme = godunov::scheme::WENOZ;
godunov::scheme mflux_scheme = godunov::scheme::UPWIND;
std::string godunov_type{"weno_z"};
std::string mflux_type{"upwind"};
const bool fluxes_are_area_weighted{false};
bool godunov_use_forces_in_trans{false};
bool m_allow_inflow_on_outflow{false};
Expand Down
1 change: 0 additions & 1 deletion amr-wind/equation_systems/icns/icns_advection.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
#include "amr-wind/equation_systems/AdvOp_Godunov.H"
#include "amr-wind/equation_systems/AdvOp_MOL.H"
#include "amr-wind/equation_systems/icns/icns.H"
#include "amr-wind/equation_systems/vof/vof_momentum_flux.H"
#include "amr-wind/core/Physics.H"

#include "AMReX_MultiFabUtil.H"
Expand Down
2 changes: 1 addition & 1 deletion amr-wind/equation_systems/temperature/temperature.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ struct Temperature : ScalarTransport
static constexpr amrex::Real default_bc_value = 0.0;

static constexpr int ndim = 1;
static constexpr bool multiply_rho = false;
static constexpr bool multiply_rho = true;
static constexpr bool has_diffusion = true;
static constexpr bool need_nph_state = true;
};
Expand Down
11 changes: 7 additions & 4 deletions amr-wind/equation_systems/vof/vof_momentum_flux.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@

namespace amr_wind::multiphase {

static void hybrid_fluxes(
static inline void hybrid_fluxes(
const FieldRepo& repo,
const int ncomp,
const amrex::Gpu::DeviceVector<int>& iconserv,
Expand All @@ -27,7 +27,8 @@ static void hybrid_fluxes(
const amrex::Real dt,
godunov::scheme mflux_scheme,
bool allow_inflow_on_outflow,
bool use_forces_in_trans)
bool use_forces_in_trans,
bool pre_multiplied_src_term = false)
{
// Get geometry
const auto& geom = repo.mesh().Geom();
Expand Down Expand Up @@ -75,8 +76,10 @@ static void hybrid_fluxes(
for (int idim = 0; idim < dof_field.num_comp(); ++idim) {
amrex::MultiFab::Multiply(
q, rho_o(lev), 0, idim, 1, fvm::Godunov::nghost_state);
amrex::MultiFab::Multiply(
fq, rho_o(lev), 0, idim, 1, fvm::Godunov::nghost_src);
if (!pre_multiplied_src_term) {
amrex::MultiFab::Multiply(
fq, rho_o(lev), 0, idim, 1, fvm::Godunov::nghost_src);
}
}

amrex::MultiFab frho(
Expand Down
16 changes: 3 additions & 13 deletions amr-wind/incflo_advance.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -297,19 +297,6 @@ void incflo::ApplyPredictor(
: amr_wind::FieldState::NPH;
icns().pre_advection_actions(fstate_preadv);

// For scalars only first
// *************************************************************************************
// if ( m_use_godunov) Compute the explicit advective terms
// R_u^(n+1/2), R_s^(n+1/2) and R_t^(n+1/2)
// if (!m_use_godunov) Compute the explicit advective terms
// R_u^n , R_s^n and R_t^n
// *************************************************************************************
// TODO: Advection computation for scalar equations have not been adjusted
// for mesh mapping
for (auto& seqn : scalar_eqns()) {
seqn->compute_advection_term(amr_wind::FieldState::Old);
}

// *************************************************************************************
// Update density first
// *************************************************************************************
Expand All @@ -323,6 +310,9 @@ void incflo::ApplyPredictor(
// updated density at `n+1/2` to be computed before other scalars use it
// when computing their source terms.
for (auto& eqn : scalar_eqns()) {
// Compute explicit advection
eqn->compute_advection_term(amr_wind::FieldState::Old);

// Compute (recompute for Godunov) the scalar forcing terms
eqn->compute_source_term(amr_wind::FieldState::NPH);

Expand Down