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

Clean up advection scheme blending #2097

Merged
merged 8 commits into from
Feb 4, 2025
Merged

Clean up advection scheme blending #2097

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
e785127
Mention upwinding fraction in the docs
ewquon Feb 3, 2025
a855498
Add note about building with homebrew openmpi
ewquon Feb 4, 2025
08e4e34
Make code formatting consistent between ERF_AdvectionSrcForMom_*.cpp
ewquon Feb 4, 2025
b7b6f08
Save upw_frac as member var (UPWIND3,UPWIND5); cleanup interp structs
ewquon Feb 4, 2025
88b9bdf
Update instantiation and calls to interp structs
ewquon Feb 4, 2025
524290d
Merge branch 'development' into upw_frac_cleanup
ewquon Feb 4, 2025
f545947
Merge branch 'development' into upw_frac_cleanup
AMLattanzi Feb 4, 2025
cad5561
Merge branch 'development' into upw_frac_cleanup
AMLattanzi Feb 4, 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
4 changes: 4 additions & 0 deletions Docs/sphinx_doc/Inputs.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -866,6 +866,10 @@ The allowed advection types for the dry and moist scalars are
"Centered_6th" and in addition,
"WENO3", "WENOZ3", "WENOMZQ3", "WENO5", "WENOZ5", "WENO7", and "WENOZ7."

To use the blended schemes, ``erf.[vartype]_[horiz|vert]_upw_frac`` for the corresponding
``erf.[vartype]_[horiz|vert]_adv_type`` should be set to a scalar between 0 and 1, where 1 is fully
upwind and 0 is fully central.

Note: if using WENO schemes, the horizontal and vertical advection types must be set to
the same string.

Expand Down
2 changes: 2 additions & 0 deletions Docs/sphinx_doc/building.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,6 +202,8 @@ Mac with CMake
~~~~~~~~~~~~~~
Tested with macOS 12.7 (Monterey) using cmake (3.27.8), open-mpi (5.0.0), and
pkg-config (0.29.2) installed with the homebrew package manager.
(Note: The homebrew openmpi library was "poured from bottle," not installed
with the ``--build-from-source`` option.)
NetCDF will be compiled from source. The instructions below should be version
agnostic.

Expand Down
48 changes: 24 additions & 24 deletions Source/Advection/ERF_AdvectionSrcForMom_ConstantDz.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,9 +54,14 @@ AdvectionSrcForMom_ConstantDz (const Box& bxx, const Box& bxy, const Box& bxz,
{
BL_PROFILE_VAR("AdvectionSrcForMom_ConstantDz", AdvectionSrcForMom_ConstantDz);

AMREX_ALWAYS_ASSERT(bxz.smallEnd(2) > 0);

auto dxInv = cellSizeInv[0], dyInv = cellSizeInv[1], dzInv = cellSizeInv[2];

AMREX_ALWAYS_ASSERT(bxz.smallEnd(2) > 0);
const bool use_terrain_fitted_coords = ( terrain_type == TerrainType::StaticFittedMesh ||
terrain_type == TerrainType::MovingFittedMesh);

AMREX_ALWAYS_ASSERT(!use_terrain_fitted_coords && (terrain_type != TerrainType::EB));

// compute mapfactor inverses
Box box2d_u(bxx); box2d_u.setRange(2,0); box2d_u.grow({3,3,0});
Expand All @@ -66,11 +71,6 @@ AdvectionSrcForMom_ConstantDz (const Box& bxx, const Box& bxy, const Box& bxz,
const Array4<Real>& mf_u_inv = mf_u_invFAB.array();
const Array4<Real>& mf_v_inv = mf_v_invFAB.array();

const bool use_terrain_fitted_coords = ( terrain_type == TerrainType::StaticFittedMesh ||
terrain_type == TerrainType::MovingFittedMesh);

AMREX_ALWAYS_ASSERT(!use_terrain_fitted_coords && (terrain_type != TerrainType::EB));

ParallelFor(box2d_u, box2d_v,
[=] AMREX_GPU_DEVICE (int i, int j, int) noexcept
{
Expand Down Expand Up @@ -146,12 +146,12 @@ AdvectionSrcForMom_ConstantDz (const Box& bxx, const Box& bxy, const Box& bxz,
} else {
if (horiz_adv_type == AdvType::Centered_2nd) {
AdvectionSrcForMomVert_N<CENTERED2>(bxx, bxy, bxz,
rho_u_rhs, rho_v_rhs, rho_w_rhs,
rho_u, rho_v, omega, u, v, w,
cellSizeInv, stretched_dz_d, mf_m,
mf_u_inv, mf_v_inv,
horiz_upw_frac, vert_upw_frac,
vert_adv_type, lo_z_face, hi_z_face);
rho_u_rhs, rho_v_rhs, rho_w_rhs,
rho_u, rho_v, omega, u, v, w,
cellSizeInv, stretched_dz_d, mf_m,
mf_u_inv, mf_v_inv,
horiz_upw_frac, vert_upw_frac,
vert_adv_type, lo_z_face, hi_z_face);
} else if (horiz_adv_type == AdvType::Upwind_3rd) {
AdvectionSrcForMomVert_N<UPWIND3>(bxx, bxy, bxz,
rho_u_rhs, rho_v_rhs, rho_w_rhs,
Expand All @@ -162,12 +162,12 @@ AdvectionSrcForMom_ConstantDz (const Box& bxx, const Box& bxy, const Box& bxz,
vert_adv_type, lo_z_face, hi_z_face);
} else if (horiz_adv_type == AdvType::Centered_4th) {
AdvectionSrcForMomVert_N<CENTERED4>(bxx, bxy, bxz,
rho_u_rhs, rho_v_rhs, rho_w_rhs,
rho_u, rho_v, omega, u, v, w,
cellSizeInv, stretched_dz_d, mf_m,
mf_u_inv, mf_v_inv,
horiz_upw_frac, vert_upw_frac,
vert_adv_type, lo_z_face, hi_z_face);
rho_u_rhs, rho_v_rhs, rho_w_rhs,
rho_u, rho_v, omega, u, v, w,
cellSizeInv, stretched_dz_d, mf_m,
mf_u_inv, mf_v_inv,
horiz_upw_frac, vert_upw_frac,
vert_adv_type, lo_z_face, hi_z_face);
} else if (horiz_adv_type == AdvType::Upwind_5th) {
AdvectionSrcForMomVert_N<UPWIND5>(bxx, bxy, bxz,
rho_u_rhs, rho_v_rhs, rho_w_rhs,
Expand All @@ -178,12 +178,12 @@ AdvectionSrcForMom_ConstantDz (const Box& bxx, const Box& bxy, const Box& bxz,
vert_adv_type, lo_z_face, hi_z_face);
} else if (horiz_adv_type == AdvType::Centered_6th) {
AdvectionSrcForMomVert_N<CENTERED6>(bxx, bxy, bxz,
rho_u_rhs, rho_v_rhs, rho_w_rhs,
rho_u, rho_v, omega, u, v, w,
cellSizeInv, stretched_dz_d, mf_m,
mf_u_inv, mf_v_inv,
horiz_upw_frac, vert_upw_frac,
vert_adv_type, lo_z_face, hi_z_face);
rho_u_rhs, rho_v_rhs, rho_w_rhs,
rho_u, rho_v, omega, u, v, w,
cellSizeInv, stretched_dz_d, mf_m,
mf_u_inv, mf_v_inv,
horiz_upw_frac, vert_upw_frac,
vert_adv_type, lo_z_face, hi_z_face);
} else {
AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
}
Expand Down
65 changes: 28 additions & 37 deletions Source/Advection/ERF_AdvectionSrcForMom_N.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,6 @@ AdvectionSrcForXMom_N (int i, int j, int k,
const amrex::Array4<const amrex::Real>& rho_w,
InterpType_H interp_u_h,
InterpType_V interp_u_v,
const amrex::Real upw_frac_h,
const amrex::Real upw_frac_v,
const amrex::Array4<const amrex::Real>& mf_u_inv,
const amrex::Array4<const amrex::Real>& mf_v_inv,
const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
Expand All @@ -43,27 +41,27 @@ AdvectionSrcForXMom_N (int i, int j, int k,
amrex::Real interp_hi(0.), interp_lo(0.);

rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h);
interp_u_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
xflux_hi = rho_u_avg_hi * interp_hi;

rho_u_avg_lo = 0.5 * (rho_u(i-1, j, k) * mf_u_inv(i-1,j,0) + rho_u(i, j, k) * mf_u_inv(i,j,0));
interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h);
interp_u_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
xflux_lo = rho_u_avg_lo * interp_lo;

rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) * mf_v_inv(i,j+1,0) + rho_v(i-1, j+1, k) * mf_v_inv(i-1,j+1,0));
interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h);
interp_u_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
yflux_hi = rho_v_avg_hi * interp_hi;

rho_v_avg_lo = 0.5 * (rho_v(i, j , k) * mf_v_inv(i,j ,0) + rho_v(i-1, j , k) * mf_v_inv(i-1,j ,0));
interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h);
interp_u_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
yflux_lo = rho_v_avg_lo * interp_lo;

rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i-1, j, k+1));
interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,upw_frac_v);
interp_u_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
zflux_hi = rho_w_avg_hi * interp_hi;

rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i-1, j, k ));
interp_u_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,upw_frac_v);
interp_u_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
zflux_lo = rho_w_avg_lo * interp_lo;

amrex::Real mfsq = 1 / (mf_u_inv(i,j,0) * mf_u_inv(i,j,0));
Expand Down Expand Up @@ -98,8 +96,6 @@ AdvectionSrcForYMom_N (int i, int j, int k,
const amrex::Array4<const amrex::Real>& rho_w,
InterpType_H interp_v_h,
InterpType_V interp_v_v,
const amrex::Real upw_frac_h,
const amrex::Real upw_frac_v,
const amrex::Array4<const amrex::Real>& mf_u_inv,
const amrex::Array4<const amrex::Real>& mf_v_inv,
const amrex::Real dxInv, const amrex::Real dyInv, const amrex::Real dzInv)
Expand All @@ -117,27 +113,27 @@ AdvectionSrcForYMom_N (int i, int j, int k,
amrex::Real interp_hi(0.), interp_lo(0.);

rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) * mf_u_inv(i+1,j ,0) + rho_u(i+1, j-1, k) * mf_u_inv(i+1,j-1,0));
interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h);
interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
xflux_hi = rho_u_avg_hi * interp_hi;

rho_u_avg_lo = 0.5 * (rho_u(i , j, k) * mf_u_inv(i ,j ,0) + rho_u(i , j-1, k) * mf_u_inv(i ,j-1,0));
interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h);
interp_v_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
xflux_lo = rho_u_avg_lo * interp_lo;

rho_v_avg_hi = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j+1, k) * mf_v_inv(i ,j+1,0));
interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h);
interp_v_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
yflux_hi = rho_v_avg_hi * interp_hi;

rho_v_avg_lo = 0.5 * (rho_v(i, j, k) * mf_v_inv(i ,j ,0) + rho_v(i, j-1, k) * mf_v_inv(i ,j-1,0));
interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h);
interp_v_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
yflux_lo = rho_v_avg_lo * interp_lo;

rho_w_avg_hi = 0.5 * (rho_w(i, j, k+1) + rho_w(i, j-1, k+1));
interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,upw_frac_v);
interp_v_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
zflux_hi = rho_w_avg_hi * interp_hi;

rho_w_avg_lo = 0.5 * (rho_w(i, j, k ) + rho_w(i, j-1, k ));
interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,rho_w_avg_lo,upw_frac_v);
interp_v_v.InterpolateInZ(i,j,k ,0,interp_lo,rho_w_avg_lo);
zflux_lo = rho_w_avg_lo * interp_lo;

amrex::Real mfsq = 1 / (mf_v_inv(i,j,0) * mf_v_inv(i,j,0));
Expand Down Expand Up @@ -176,8 +172,6 @@ AdvectionSrcForZMom_N (int i, int j, int k,
InterpType_H interp_w_h,
InterpType_V interp_w_v,
WallInterpType interp_w_wall,
const amrex::Real upw_frac_h,
const amrex::Real upw_frac_v,
const amrex::Array4<const amrex::Real>& mf_m,
const amrex::Array4<const amrex::Real>& mf_u_inv,
const amrex::Array4<const amrex::Real>& mf_v_inv,
Expand All @@ -199,19 +193,19 @@ AdvectionSrcForZMom_N (int i, int j, int k,
amrex::Real interp_hi(0.), interp_lo(0.);

rho_u_avg_hi = 0.5 * (rho_u(i+1, j, k) + rho_u(i+1, j, k-1)) * mf_u_inv(i+1,j ,0);
interp_w_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h);
interp_w_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi);
xflux_hi = rho_u_avg_hi * interp_hi;

rho_u_avg_lo = 0.5 * (rho_u(i , j, k) + rho_u(i , j, k-1)) * mf_u_inv(i ,j ,0);
interp_w_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h);
interp_w_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo);
xflux_lo = rho_u_avg_lo * interp_lo;

rho_v_avg_hi = 0.5 * (rho_v(i, j+1, k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0);
interp_w_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi,upw_frac_h);
interp_w_h.InterpolateInY(i,j+1,k,0,interp_hi,rho_v_avg_hi);
yflux_hi = rho_v_avg_hi * interp_hi;

rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0);
interp_w_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h);
interp_w_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo);
yflux_lo = rho_v_avg_lo * interp_lo;

// int l_spatial_order_hi = std::min(std::min(vert_spatial_order, 2*(hi_z_face-k)), 2*(k+1));
Expand All @@ -225,15 +219,15 @@ AdvectionSrcForZMom_N (int i, int j, int k,
rho_w_avg_hi = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k+1));
if (k == hi_z_face-1)
{
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,0,AdvType::Centered_2nd);
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_2nd);
} else if (k == hi_z_face-2 || k == lo_z_face+1) {
if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,0,AdvType::Centered_4th);
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,AdvType::Centered_4th);
} else {
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,upw_frac_v,vert_adv_type);
interp_w_wall.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,vert_adv_type);
}
} else {
interp_w_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi,upw_frac_v);
interp_w_v.InterpolateInZ(i,j,k+1,0,interp_hi,rho_w_avg_hi);
}
zflux_hi = rho_w_avg_hi * interp_hi;
}
Expand All @@ -248,15 +242,15 @@ AdvectionSrcForZMom_N (int i, int j, int k,
} else {
rho_w_avg_lo = 0.5 * (rho_w(i,j,k) + rho_w(i,j,k-1));
if (k == lo_z_face+1) {
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,0,AdvType::Centered_2nd);
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_2nd);
} else if (k == lo_z_face+2 || k == hi_z_face-1) {
if (vert_adv_type != AdvType::Centered_2nd && vert_adv_type != AdvType::Upwind_3rd) {
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,0,AdvType::Centered_4th);
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,AdvType::Centered_4th);
} else {
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,upw_frac_v,vert_adv_type);
interp_w_wall.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,vert_adv_type);
}
} else {
interp_w_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo,upw_frac_v);
interp_w_v.InterpolateInZ(i,j,k,0,interp_lo,rho_w_avg_lo);
}
zflux_lo = rho_w_avg_lo * interp_lo;
}
Expand Down Expand Up @@ -296,10 +290,10 @@ AdvectionSrcForMomWrapper_N (const amrex::Box& bxx, const amrex::Box& bxy, const
const int lo_z_face, const int hi_z_face)
{
// Instantiate the appropriate structs
InterpType_H interp_u_h(u); InterpType_V interp_u_v(u); // X-MOM
InterpType_H interp_v_h(v); InterpType_V interp_v_v(v); // Y-MOM
InterpType_H interp_w_h(w); InterpType_V interp_w_v(w); // Z-MOM
WallInterpType interp_w_wall(w); // Z-MOM @ wall
InterpType_H interp_u_h(u, upw_frac_h); InterpType_V interp_u_v(u, upw_frac_v); // X-MOM
InterpType_H interp_v_h(v, upw_frac_h); InterpType_V interp_v_v(v, upw_frac_v); // Y-MOM
InterpType_H interp_w_h(w, upw_frac_h); InterpType_V interp_w_v(w, upw_frac_v); // Z-MOM
WallInterpType interp_w_wall(w, upw_frac_v); // Z-MOM @ wall

auto dxInv = cellSizeInv[0];
auto dyInv = cellSizeInv[1];
Expand All @@ -312,7 +306,6 @@ AdvectionSrcForMomWrapper_N (const amrex::Box& bxx, const amrex::Box& bxy, const
auto dzInv = 1.0/dz_ptr[k];
rho_u_rhs(i, j, k) = -AdvectionSrcForXMom_N(i, j, k, rho_u, rho_v, rho_w,
interp_u_h, interp_u_v,
upw_frac_h, upw_frac_v,
mf_u_inv, mf_v_inv,
dxInv, dyInv, dzInv);
});
Expand All @@ -323,7 +316,6 @@ AdvectionSrcForMomWrapper_N (const amrex::Box& bxx, const amrex::Box& bxy, const
auto dzInv = 1.0/dz_ptr[k];
rho_v_rhs(i, j, k) = -AdvectionSrcForYMom_N(i, j, k, rho_u, rho_v, rho_w,
interp_v_h, interp_v_v,
upw_frac_h, upw_frac_v,
mf_u_inv, mf_v_inv,
dxInv, dyInv, dzInv);
});
Expand All @@ -334,7 +326,6 @@ AdvectionSrcForMomWrapper_N (const amrex::Box& bxx, const amrex::Box& bxy, const
auto dzInv = 2.0/(dz_ptr[k] + dz_ptr[k-1]);
rho_w_rhs(i, j, k) = -AdvectionSrcForZMom_N(i, j, k, rho_u, rho_v, rho_w, w,
interp_w_h, interp_w_v, interp_w_wall,
upw_frac_h, upw_frac_v,
mf_m, mf_u_inv, mf_v_inv,
vert_adv_type, lo_z_face, hi_z_face,
dxInv, dyInv, dzInv);
Expand Down
Loading
Loading