Skip to content

Commit

Permalink
Make numdiff match WRF docs. (#1929)
Browse files Browse the repository at this point in the history 
* Make numdiff match WRF docs.

* comment unused vars.

* one mapfac multiply.

* update ctest.
  • Loading branch information
@AMLattanzi
AMLattanzi authored Nov 5, 2024
1 parent 9428c70 commit d55d717
Show file tree
Hide file tree
Showing 10 changed files with 247 additions and 153 deletions.
59 changes: 41 additions & 18 deletions Source/SourceTerms/ERF_NumericalDiffusion.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,23 +5,46 @@
#include <ERF_DataStruct.H>
#include <AMReX_MultiFab.H>

void NumericalDiffusion (const amrex::Box& bx,
const int start_comp,
const int num_comp,
const amrex::Real dt,
const amrex::Real num_diff_coeff,
const amrex::Array4<const amrex::Real>& data,
const amrex::Array4< amrex::Real>& rhs,
const amrex::Array4<const amrex::Real>& mf_x,
const amrex::Array4<const amrex::Real>& mf_y,
const bool avg_mf_x_y,
const bool avg_mf_y_x);

void NumericalDiffusionVert (const amrex::Box& bx,
const int start_comp,
const int num_comp,
const amrex::Real dt,
const amrex::Real num_diff_coeff,
const amrex::Array4<const amrex::Real>& data,
const amrex::Array4< amrex::Real>& rhs);
AMREX_GPU_DEVICE
AMREX_FORCE_INLINE
amrex::Real
calc_fifth_order_deriv (const amrex::Real& dnp2,
const amrex::Real& dnp1,
const amrex::Real& dn,
const amrex::Real& dnm1,
const amrex::Real& dnm2,
const amrex::Real& dnm3)
{
amrex::Real interp = 10. * (dn - dnm1)
- 5. * (dnp1 - dnm2)
+ (dnp2 - dnm3);
return interp;
}

void NumericalDiffusion_Scal (const amrex::Box& bx,
const int start_comp,
const int num_comp,
const amrex::Real dt,
const amrex::Real num_diff_coeff,
const amrex::Array4<const amrex::Real>& prim_data,
const amrex::Array4<const amrex::Real>& cell_data,
const amrex::Array4< amrex::Real>& rhs,
const amrex::Array4<const amrex::Real>& mf);

void NumericalDiffusion_Xmom (const amrex::Box& bx,
const amrex::Real dt,
const amrex::Real num_diff_coeff,
const amrex::Array4<const amrex::Real>& prim_data,
const amrex::Array4<const amrex::Real>& cell_data,
const amrex::Array4< amrex::Real>& rhs,
const amrex::Array4<const amrex::Real>& mf);

void NumericalDiffusion_Ymom (const amrex::Box& bx,
const amrex::Real dt,
const amrex::Real num_diff_coeff,
const amrex::Array4<const amrex::Real>& prim_data,
const amrex::Array4<const amrex::Real>& cell_data,
const amrex::Array4< amrex::Real>& rhs,
const amrex::Array4<const amrex::Real>& mf);
#endif
228 changes: 153 additions & 75 deletions Source/SourceTerms/ERF_NumericalDiffusion.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,112 +9,190 @@ using namespace amrex;
* @param[in] start_comp staring component index
* @param[in] num_comp number of total components
* @param[in] num_diff_coeff
* @param[in] data variables used to compute RHS
* @param[in] prim_data primitive variables
* @param[in] cell_data cell center variables
* @param[out] rhs store the right hand side
* @param[in] mf_x map factor at x-face
* @param[in] mf_y map factor at y-face
* @param[in] avg_mf_x_y flag to average map factor x in y-dir
* @param[in] avg_mf_y_x flag to average map factor y in x-dir
* @param[in] mf map factor
*/
void
NumericalDiffusion (const Box& bx,
const int start_comp,
const int num_comp,
const Real dt,
const Real num_diff_coeff,
const Array4<const Real>& data,
const Array4< Real>& rhs,
const Array4<const Real>& mf_x,
const Array4<const Real>& mf_y,
const bool avg_mf_x_y,
const bool avg_mf_y_x)
NumericalDiffusion_Scal (const Box& bx,
const int start_comp,
const int num_comp,
const Real dt,
const Real num_diff_coeff,
const Array4<const Real>& prim_data,
const Array4<const Real>& cell_data,
const Array4< Real>& rhs,
const Array4<const Real>& mf_arr)
{
BL_PROFILE_VAR("NumericalDiffusion()",NumericalDiffusion);

// Average map factors to correct locations
Box planebx(bx); planebx.setSmall(2,0); planebx.setBig(2,0);
FArrayBox mf_x_bar; mf_x_bar.resize(planebx,1,The_Async_Arena());
FArrayBox mf_y_bar; mf_y_bar.resize(planebx,1,The_Async_Arena());
const Array4<Real>& mfx_arr = mf_x_bar.array();
const Array4<Real>& mfy_arr = mf_y_bar.array();
ParallelFor(planebx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
BL_PROFILE_VAR("NumericalDiffusion_Scal()",NumericalDiffusion_Scal);

// Capture diffusion coeff
Real coeff6 = num_diff_coeff / (2.0 * dt);

// Compute 5th order derivative and augment RHS
ParallelFor(bx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int m) noexcept
{
if (avg_mf_x_y) {
mfx_arr(i,j,k) = 0.5 * ( mf_x(i,j-1,k) + mf_x(i,j,k) );
} else {
mfx_arr(i,j,k) = mf_x(i,j,k);
}
if (avg_mf_y_x) {
mfy_arr(i,j,k) = 0.5 * ( mf_y(i-1,j,k) + mf_y(i,j,k) );
} else {
mfy_arr(i,j,k) = mf_y(i,j,k);
}
int n = start_comp + m; // conserved index
int nm1 = (n==0) ? 0 : n - 1; // prim index
Real rho_x_lo = (n==0) ? 1.0 : 0.5 * ( cell_data(i-1,j,k,Rho_comp) + cell_data(i ,j,k,Rho_comp) );
Real xflux_lo = rho_x_lo * calc_fifth_order_deriv(prim_data(i+2,j,k,nm1), prim_data(i+1,j,k,nm1),
prim_data(i ,j,k,nm1), prim_data(i-1,j,k,nm1),
prim_data(i-2,j,k,nm1), prim_data(i-3,j,k,nm1));
if ( (xflux_lo * (prim_data(i,j,k,nm1) - prim_data(i-1,j,k,nm1)) ) < 0.) xflux_lo = 0.;


Real rho_x_hi = (n==0) ? 1.0 : 0.5 * ( cell_data(i+1,j,k,Rho_comp) + cell_data(i ,j,k,Rho_comp) );
Real xflux_hi = rho_x_hi * calc_fifth_order_deriv(prim_data(i+3,j,k,nm1), prim_data(i+2,j,k,nm1),
prim_data(i+1,j,k,nm1), prim_data(i ,j,k,nm1),
prim_data(i-1,j,k,nm1), prim_data(i-2,j,k,nm1));
if ( (xflux_hi * (prim_data(i+1,j,k,nm1) - prim_data(i,j,k,nm1)) ) < 0.) xflux_hi = 0.;


Real rho_y_lo = (n==0) ? 1.0 : 0.5 * ( cell_data(i,j-1,k,Rho_comp) + cell_data(i,j ,k,Rho_comp) );
Real yflux_lo = rho_y_lo * calc_fifth_order_deriv(prim_data(i,j+2,k,nm1), prim_data(i,j+1,k,nm1),
prim_data(i,j ,k,nm1), prim_data(i,j-1,k,nm1),
prim_data(i,j-2,k,nm1), prim_data(i,j-3,k,nm1));
if ( (yflux_lo * (prim_data(i,j,k,nm1) - prim_data(i,j-1,k,nm1)) ) < 0.) yflux_lo = 0.;


Real rho_y_hi = (n==0) ? 1.0 : 0.5 * ( cell_data(i,j+1,k,Rho_comp) + cell_data(i,j ,k,Rho_comp) );
Real yflux_hi = rho_y_hi * calc_fifth_order_deriv(prim_data(i,j+3,k,nm1), prim_data(i,j+2,k,nm1),
prim_data(i,j+1,k,nm1), prim_data(i,j ,k,nm1),
prim_data(i,j-1,k,nm1), prim_data(i,j-2,k,nm1));
if ( (yflux_hi * (prim_data(i,j+1,k,nm1) - prim_data(i,j,k,nm1)) ) < 0.) yflux_hi = 0.;


rhs(i,j,k,n) += coeff6 * mf_arr(i,j,0) * ( (xflux_hi - xflux_lo)
+ (yflux_hi - yflux_lo) );
});
}

/**
* Function to compute 6th order numerical diffusion RHS.
*
* @param[in] bx box to loop over
* @param[in] start_comp staring component index
* @param[in] num_comp number of total components
* @param[in] num_diff_coeff
* @param[in] prim_data primitive variables
* @param[in] cell_data cell center variables
* @param[out] rhs store the right hand side
* @param[in] mf map factor
*/
void
NumericalDiffusion_Xmom (const Box& bx,
const Real dt,
const Real num_diff_coeff,
const Array4<const Real>& prim_data,
const Array4<const Real>& cell_data,
const Array4< Real>& rhs,
const Array4<const Real>& mf_arr)
{
BL_PROFILE_VAR("NumericalDiffusion_Xmom()",NumericalDiffusion_Xmom);

// Capture diffusion coeff
Real coeff6 = num_diff_coeff / (2.0 * dt);

// Compute 5th order derivative and augment RHS
ParallelFor(bx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int m) noexcept
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
int n = start_comp + m;
Real xflux_lo = 10. * (data(i ,j,k,n) - data(i-1,j,k,n))
- 5. * (data(i+1,j,k,n) - data(i-2,j,k,n))
+ (data(i+2,j,k,n) - data(i-3,j,k,n));
if ( (xflux_lo * (data(i,j,k,n) - data(i-1,j,k,n)) ) < 0.) xflux_lo = 0.;
Real xflux_hi = 10. * (data(i+1,j,k,n) - data(i ,j,k,n))
- 5. * (data(i+2,j,k,n) - data(i-1,j,k,n))
+ (data(i+3,j,k,n) - data(i-2,j,k,n));
if ( (xflux_hi * (data(i+1,j,k,n) - data(i,j,k,n)) ) < 0.) xflux_hi = 0.;
Real yflux_lo = 10. * (data(i,j ,k,n) - data(i,j-1,k,n))
- 5. * (data(i,j+1,k,n) - data(i,j-2,k,n))
+ (data(i,j+2,k,n) - data(i,j-3,k,n));
if ( (yflux_lo * (data(i,j,k,n) - data(i,j-1,k,n)) ) < 0.) yflux_lo = 0.;
Real yflux_hi = 10. * (data(i,j+1,k,n) - data(i,j ,k,n))
- 5. * (data(i,j+2,k,n) - data(i,j-1,k,n))
+ (data(i,j+3,k,n) - data(i,j-2,k,n));
if ( (yflux_hi * (data(i,j+1,k,n) - data(i,j,k,n)) ) < 0.) yflux_hi = 0.;
rhs(i,j,k,n) += coeff6 * ( (xflux_hi - xflux_lo) * mfx_arr(i,j,0)
+ (yflux_hi - yflux_lo) * mfy_arr(i,j,0) );
Real rho_x_lo = cell_data(i-1,j,k,Rho_comp);
Real xflux_lo = rho_x_lo * calc_fifth_order_deriv(prim_data(i+2,j,k), prim_data(i+1,j,k),
prim_data(i ,j,k), prim_data(i-1,j,k),
prim_data(i-2,j,k), prim_data(i-3,j,k));
if ( (xflux_lo * (prim_data(i,j,k) - prim_data(i-1,j,k)) ) < 0.) xflux_lo = 0.;


Real rho_x_hi = cell_data(i ,j,k,Rho_comp);
Real xflux_hi = rho_x_hi * calc_fifth_order_deriv(prim_data(i+3,j,k), prim_data(i+2,j,k),
prim_data(i+1,j,k), prim_data(i ,j,k),
prim_data(i-1,j,k), prim_data(i-2,j,k));
if ( (xflux_hi * (prim_data(i+1,j,k) - prim_data(i,j,k)) ) < 0.) xflux_hi = 0.;


Real rho_y_lo = 0.25 * ( cell_data(i ,j ,k,Rho_comp) + cell_data(i-1,j ,k,Rho_comp)
+ cell_data(i ,j-1,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
Real yflux_lo = rho_y_lo * calc_fifth_order_deriv(prim_data(i,j+2,k), prim_data(i,j+1,k),
prim_data(i,j ,k), prim_data(i,j-1,k),
prim_data(i,j-2,k), prim_data(i,j-3,k));
if ( (yflux_lo * (prim_data(i,j,k) - prim_data(i,j-1,k)) ) < 0.) yflux_lo = 0.;


Real rho_y_hi = 0.25 * ( cell_data(i ,j ,k,Rho_comp) + cell_data(i-1,j ,k,Rho_comp)
+ cell_data(i ,j+1,k,Rho_comp) + cell_data(i-1,j+1,k,Rho_comp) );
Real yflux_hi = rho_y_hi * calc_fifth_order_deriv(prim_data(i,j+3,k), prim_data(i,j+2,k),
prim_data(i,j+1,k), prim_data(i,j ,k),
prim_data(i,j-1,k), prim_data(i,j-2,k));
if ( (yflux_hi * (prim_data(i,j+1,k) - prim_data(i,j,k)) ) < 0.) yflux_hi = 0.;


rhs(i,j,k) += coeff6 * mf_arr(i,j,0) * ( (xflux_hi - xflux_lo)
+ (yflux_hi - yflux_lo) );
});
}


/**
* Function to compute 6th order numerical diffusion RHS.
*
* @param[in] bx box to loop over
* @param[in] start_comp staring component index
* @param[in] num_comp number of total components
* @param[in] num_diff_coeff
* @param[in] data variables used to compute RHS
* @param[in] prim_data primitive variables
* @param[in] cell_data cell center variables
* @param[out] rhs store the right hand side
* @param[in] mf map factor
*/
void
NumericalDiffusionVert (const Box& bx,
const int start_comp,
const int num_comp,
const Real dt,
const Real num_diff_coeff,
const Array4<const Real>& data,
const Array4< Real>& rhs)
NumericalDiffusion_Ymom (const Box& bx,
const Real dt,
const Real num_diff_coeff,
const Array4<const Real>& prim_data,
const Array4<const Real>& cell_data,
const Array4< Real>& rhs,
const Array4<const Real>& mf_arr)
{
BL_PROFILE_VAR("NumericalDiffusionVert()",NumericalDiffusionVert);
BL_PROFILE_VAR("NumericalDiffusion_Ymom()",NumericalDiffusion_Ymom);

// Capture diffusion coeff
Real coeff6 = num_diff_coeff / (2.0 * dt);

// Compute 5th order derivative and augment RHS
ParallelFor(bx, num_comp, [=] AMREX_GPU_DEVICE (int i, int j, int k, int m) noexcept
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
int n = start_comp + m;
Real zflux_lo = 10. * (data(i,j,k ,n) - data(i,j,k-1,n))
- 5. * (data(i,j,k+1,n) - data(i,j,k-2,n))
+ (data(i,j,k+2,n) - data(i,j,k-3,n));
if ( (zflux_lo * (data(i,j,k,n) - data(i,j,k-1,n)) ) < 0.) zflux_lo = 0.;
Real zflux_hi = 10. * (data(i,j,k+1,n) - data(i,j,k ,n))
- 5. * (data(i,j,k+2,n) - data(i,j,k-1,n))
+ (data(i,j,k+3,n) - data(i,j,k-2,n));
if ( (zflux_hi * (data(i,j,k+1,n) - data(i,j,k,n)) ) < 0.) zflux_hi = 0.;
rhs(i,j,k,n) += coeff6 * ( (zflux_hi - zflux_lo) );
Real rho_x_lo = 0.25 * ( cell_data(i ,j ,k,Rho_comp) + cell_data(i ,j-1,k,Rho_comp)
+ cell_data(i-1,j ,k,Rho_comp) + cell_data(i-1,j-1,k,Rho_comp) );
Real xflux_lo = rho_x_lo * calc_fifth_order_deriv(prim_data(i+2,j,k), prim_data(i+1,j,k),
prim_data(i ,j,k), prim_data(i-1,j,k),
prim_data(i-2,j,k), prim_data(i-3,j,k));
if ( (xflux_lo * (prim_data(i,j,k) - prim_data(i-1,j,k)) ) < 0.) xflux_lo = 0.;


Real rho_x_hi = 0.25 * ( cell_data(i ,j ,k,Rho_comp) + cell_data(i ,j-1,k,Rho_comp)
+ cell_data(i+1,j ,k,Rho_comp) + cell_data(i+1,j-1,k,Rho_comp) );
Real xflux_hi = rho_x_hi * calc_fifth_order_deriv(prim_data(i+3,j,k), prim_data(i+2,j,k),
prim_data(i+1,j,k), prim_data(i ,j,k),
prim_data(i-1,j,k), prim_data(i-2,j,k));
if ( (xflux_hi * (prim_data(i+1,j,k) - prim_data(i,j,k)) ) < 0.) xflux_hi = 0.;


Real rho_y_lo = cell_data(i,j-1,k,Rho_comp);
Real yflux_lo = rho_y_lo * calc_fifth_order_deriv(prim_data(i,j+2,k), prim_data(i,j+1,k),
prim_data(i,j ,k), prim_data(i,j-1,k),
prim_data(i,j-2,k), prim_data(i,j-3,k));
if ( (yflux_lo * (prim_data(i,j,k) - prim_data(i,j-1,k)) ) < 0.) yflux_lo = 0.;


Real rho_y_hi = cell_data(i,j ,k,Rho_comp);
Real yflux_hi = rho_y_hi * calc_fifth_order_deriv(prim_data(i,j+3,k), prim_data(i,j+2,k),
prim_data(i,j+1,k), prim_data(i,j ,k),
prim_data(i,j-1,k), prim_data(i,j-2,k));
if ( (yflux_hi * (prim_data(i,j ,k) - prim_data(i,j,k)) ) < 0.) yflux_hi = 0.;


rhs(i,j,k) += coeff6 * mf_arr(i,j,0) * ( (xflux_hi - xflux_lo)
+ (yflux_hi - yflux_lo) );
});
}
1 change: 1 addition & 0 deletions Source/SourceTerms/ERF_Src_headers.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ void make_sources (int level, int nrk,
const SolverChoice& solverChoice,
std::unique_ptr<amrex::MultiFab>& mapfac_u,
std::unique_ptr<amrex::MultiFab>& mapfac_v,
std::unique_ptr<amrex::MultiFab>& mapfac_m,
const amrex::Real* dptr_rhotheta_src,
const amrex::Real* dptr_rhoqt_src,
const amrex::Real* dptr_wbar_sub,
Expand Down
Loading

0 comments on commit d55d717

Please sign in to comment.