Skip to content

Commit

Permalink
do limiting on the interpolation of scalars.
Browse files Browse the repository at this point in the history
  • Loading branch information
@AMLattanzi
AMLattanzi committed Jul 25, 2024
1 parent bd9eda8 commit 4ff9770
Show file tree
Hide file tree
Showing 4 changed files with 47 additions and 87 deletions.
30 changes: 21 additions & 9 deletions Source/Advection/AdvectionSrcForScalars.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,10 @@ AdvectionSrcForScalarsWrapper (const amrex::Box& bx,
const amrex::Array4<const amrex::Real>& avg_ymom,
const amrex::Array4<const amrex::Real>& avg_zmom,
const amrex::Real horiz_upw_frac,
const amrex::Real vert_upw_frac)
const amrex::Real vert_upw_frac,
const bool& use_mono_adv,
amrex::Real* max_s_ptr,
amrex::Real* min_s_ptr)
{
// Instantiate structs for vert/horiz interp
InterpType_H interp_prim_h(cell_prim);
Expand All @@ -32,6 +35,7 @@ AdvectionSrcForScalarsWrapper (const amrex::Box& bx,
amrex::Real interpx(0.);

interp_prim_h.InterpolateInX(i,j,k,prim_index,interpx,avg_xmom(i ,j ,k ),horiz_upw_frac);
if (use_mono_adv) interpx = amrex::min(amrex::max(interpx,min_s_ptr[prim_index]),max_s_ptr[prim_index]);
(flx_arr[0])(i,j,k,cons_index) = avg_xmom(i,j,k) * interpx;
});
amrex::ParallelFor(ybx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -41,7 +45,7 @@ AdvectionSrcForScalarsWrapper (const amrex::Box& bx,

amrex::Real interpy(0.);
interp_prim_h.InterpolateInY(i,j,k,prim_index,interpy,avg_ymom(i ,j ,k ),horiz_upw_frac);

if (use_mono_adv) interpy = amrex::min(amrex::max(interpy,min_s_ptr[prim_index]),max_s_ptr[prim_index]);
(flx_arr[1])(i,j,k,cons_index) = avg_ymom(i,j,k) * interpy;
});
amrex::ParallelFor(zbx, ncomp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
Expand All @@ -52,7 +56,7 @@ AdvectionSrcForScalarsWrapper (const amrex::Box& bx,
amrex::Real interpz(0.);

interp_prim_v.InterpolateInZ(i,j,k,prim_index,interpz,avg_zmom(i ,j ,k ),vert_upw_frac);

if (use_mono_adv) interpz = amrex::min(amrex::max(interpz,min_s_ptr[prim_index]),max_s_ptr[prim_index]);
(flx_arr[2])(i,j,k,cons_index) = avg_zmom(i,j,k) * interpz;
});
}
Expand All @@ -71,38 +75,46 @@ AdvectionSrcForScalarsVert (const amrex::Box& bx,
const amrex::Array4<const amrex::Real>& avg_zmom,
const amrex::Real horiz_upw_frac,
const amrex::Real vert_upw_frac,
const AdvType vert_adv_type)
const AdvType vert_adv_type,
const bool& use_mono_adv,
amrex::Real* max_s_ptr,
amrex::Real* min_s_ptr)
{
switch(vert_adv_type) {
case AdvType::Centered_2nd:
AdvectionSrcForScalarsWrapper<InterpType_H,CENTERED2>(bx, ncomp, icomp,
flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Upwind_3rd:
AdvectionSrcForScalarsWrapper<InterpType_H,UPWIND3>(bx, ncomp, icomp,
flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Centered_4th:
AdvectionSrcForScalarsWrapper<InterpType_H,CENTERED4>(bx, ncomp, icomp,
flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Upwind_5th:
AdvectionSrcForScalarsWrapper<InterpType_H,UPWIND5>(bx, ncomp, icomp,
flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Centered_6th:
AdvectionSrcForScalarsWrapper<InterpType_H,CENTERED6>(bx, ncomp, icomp,
flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
default:
AMREX_ASSERT_WITH_MESSAGE(false, "Unknown vertical advection scheme!");
Expand Down
92 changes: 20 additions & 72 deletions Source/Advection/AdvectionSrcForState.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -202,120 +202,68 @@ AdvectionSrcForScalars (const Real& dt,
case AdvType::Centered_2nd:
AdvectionSrcForScalarsVert<CENTERED2>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac, vert_adv_type);
horiz_upw_frac, vert_upw_frac, vert_adv_type,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Upwind_3rd:
AdvectionSrcForScalarsVert<UPWIND3>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac, vert_adv_type);
horiz_upw_frac, vert_upw_frac, vert_adv_type,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Centered_4th:
AdvectionSrcForScalarsVert<CENTERED4>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac, vert_adv_type);
horiz_upw_frac, vert_upw_frac, vert_adv_type,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Upwind_5th:
AdvectionSrcForScalarsVert<UPWIND5>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac, vert_adv_type);
horiz_upw_frac, vert_upw_frac, vert_adv_type,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Centered_6th:
AdvectionSrcForScalarsVert<CENTERED6>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac, vert_adv_type);
horiz_upw_frac, vert_upw_frac, vert_adv_type,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Weno_3:
AdvectionSrcForScalarsWrapper<WENO3,WENO3>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Weno_5:
AdvectionSrcForScalarsWrapper<WENO5,WENO5>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Weno_3Z:
AdvectionSrcForScalarsWrapper<WENO_Z3,WENO_Z3>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Weno_3MZQ:
AdvectionSrcForScalarsWrapper<WENO_MZQ3,WENO_MZQ3>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
case AdvType::Weno_5Z:
AdvectionSrcForScalarsWrapper<WENO_Z5,WENO_Z5>(bx, ncomp, icomp, flx_arr, cell_prim,
avg_xmom, avg_ymom, avg_zmom,
horiz_upw_frac, vert_upw_frac);
horiz_upw_frac, vert_upw_frac,
use_mono_adv, max_s_ptr, min_s_ptr);
break;
default:
AMREX_ASSERT_WITH_MESSAGE(false, "Unknown advection scheme!");
}
}

// Monotonicity preserving order reduction for SLOW SCALARS (0-th upwind)
if (use_mono_adv) {
ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
const int cons_index = icomp + n;
const int prim_index = cons_index - 1;

Real max_val = max_s_ptr[cons_index];
Real min_val = min_s_ptr[cons_index];

Real invdetJ = (detJ(i,j,k) > 0.) ? 1. / detJ(i,j,k) : 1.;
Real mfsq = mf_m(i,j,0) * mf_m(i,j,0);

Real RHS = - invdetJ * mfsq * (
( (flx_arr[0])(i+1,j ,k ,cons_index) - (flx_arr[0])(i,j,k,cons_index) ) * dxInv +
( (flx_arr[1])(i ,j+1,k ,cons_index) - (flx_arr[1])(i,j,k,cons_index) ) * dyInv +
( (flx_arr[2])(i ,j ,k+1,cons_index) - (flx_arr[2])(i,j,k,cons_index) ) * dzInv );

// NOTE: This forward prediction uses the `cur_cons` as opposed to `old_cons` since
// source terms may cause increase/decrease in a variable each RK stage. We
// want to ensure the advection operator does not induce over/under-shoot
// from the current state. If the `old_cons` is used and significant forcing is
// present, we could trip an order reduction just due to the source terms.
Real tmp_upd = cur_cons(i,j,k,cons_index) + RHS*dt;
if (tmp_upd<min_val || tmp_upd>max_val) {
// HI
if (avg_xmom(i+1,j,k)>0.0) {
(flx_arr[0])(i+1,j,k,cons_index) = avg_xmom(i+1,j,k) * cell_prim(i ,j,k,prim_index);
} else {
(flx_arr[0])(i+1,j,k,cons_index) = avg_xmom(i+1,j,k) * cell_prim(i+1,j,k,prim_index);
}
if (avg_ymom(i,j+1,k)>0.0) {
(flx_arr[1])(i,j+1,k,cons_index) = avg_ymom(i,j+1,k) * cell_prim(i,j ,k,prim_index);
} else {
(flx_arr[1])(i,j+1,k,cons_index) = avg_ymom(i,j+1,k) * cell_prim(i,j+1,k,prim_index);
}
if (avg_zmom(i,j,k+1)>0.0) {
(flx_arr[2])(i,j,k+1,cons_index) = avg_zmom(i,j,k+1) * cell_prim(i,j,k ,prim_index);
} else {
(flx_arr[2])(i,j,k+1,cons_index) = avg_zmom(i,j,k+1) * cell_prim(i,j,k+1,prim_index);
}

// LO
if (avg_xmom(i,j,k)>0.0) {
(flx_arr[0])(i,j,k,cons_index) = avg_xmom(i,j,k) * cell_prim(i-1,j,k,prim_index);
} else {
(flx_arr[0])(i,j,k,cons_index) = avg_xmom(i,j,k) * cell_prim(i ,j,k,prim_index);
}
if (avg_ymom(i,j,k)>0.0) {
(flx_arr[1])(i,j,k,cons_index) = avg_ymom(i,j,k) * cell_prim(i,j-1,k,prim_index);
} else {
(flx_arr[1])(i,j,k,cons_index) = avg_ymom(i,j,k) * cell_prim(i,j ,k,prim_index);
}
if (avg_zmom(i,j,k)>0.0) {
(flx_arr[2])(i,j,k,cons_index) = avg_zmom(i,j,k) * cell_prim(i,j,k-1,prim_index);
} else {
(flx_arr[2])(i,j,k,cons_index) = avg_zmom(i,j,k) * cell_prim(i,j,k ,prim_index);
}
}
});
}

ParallelFor(bx, ncomp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
{
Real invdetJ = (detJ(i,j,k) > 0.) ? 1. / detJ(i,j,k) : 1.;
Expand Down
6 changes: 3 additions & 3 deletions Source/TimeIntegration/ERF_slow_rhs_post.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -168,15 +168,15 @@ void erf_slow_rhs_post (int level, int finest_level,
// *****************************************************************************
// Monotonic advection for scalars
// *****************************************************************************
int nvar = S_new[IntVars::cons].nComp();
int nvar = S_prim.nComp();
Vector<Real> max_scal(nvar, 1.0e34); Gpu::DeviceVector<Real> max_scal_d(nvar);
Vector<Real> min_scal(nvar,-1.0e34); Gpu::DeviceVector<Real> min_scal_d(nvar);
if (l_use_mono_adv) {
auto const& ma_s_arr = S_new[IntVars::cons].const_arrays();
auto const& ma_s_arr = S_prim.const_arrays();
for (int ivar(RhoKE_comp); ivar<nvar; ++ivar) {
GpuTuple<Real,Real> mm = ParReduce(TypeList<ReduceOpMax,ReduceOpMin>{},
TypeList<Real, Real>{},
S_new[IntVars::cons], IntVect(0),
S_prim, S_prim.nGrowVect(),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> GpuTuple<Real,Real>
{
Expand Down
6 changes: 3 additions & 3 deletions Source/TimeIntegration/ERF_slow_rhs_pre.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,15 +203,15 @@ void erf_slow_rhs_pre (int level, int finest_level,
// *****************************************************************************
// Monotonic advection for scalars
// *****************************************************************************
int nvar = S_data[IntVars::cons].nComp();
int nvar = S_prim.nComp();
Vector<Real> max_scal(nvar, 1.0e34); Gpu::DeviceVector<Real> max_scal_d(nvar);
Vector<Real> min_scal(nvar,-1.0e34); Gpu::DeviceVector<Real> min_scal_d(nvar);
if (l_use_mono_adv) {
auto const& ma_s_arr = S_data[IntVars::cons].const_arrays();
auto const& ma_s_arr = S_prim.const_arrays();
for (int ivar(RhoTheta_comp); ivar<RhoKE_comp; ++ivar) {
GpuTuple<Real,Real> mm = ParReduce(TypeList<ReduceOpMax,ReduceOpMin>{},
TypeList<Real, Real>{},
S_data[IntVars::cons], IntVect(0),
S_prim, S_prim.nGrowVect(),
[=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept
-> GpuTuple<Real,Real>
{
Expand Down

0 comments on commit 4ff9770

Please sign in to comment.