From ea56a63767a27f770bf775491d762d9a5e9c2517 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Tue, 5 Nov 2024 12:24:37 -0800 Subject: [PATCH] Make numdiff match WRF docs. --- Source/SourceTerms/ERF_NumericalDiffusion.H | 59 +++-- Source/SourceTerms/ERF_NumericalDiffusion.cpp | 228 ++++++++++++------ Source/SourceTerms/ERF_Src_headers.H | 1 + Source/SourceTerms/ERF_make_mom_sources.cpp | 23 +- Source/SourceTerms/ERF_make_sources.cpp | 37 +-- Source/TimeIntegration/ERF_TI_slow_rhs_fun.H | 4 +- 6 files changed, 229 insertions(+), 123 deletions(-) diff --git a/Source/SourceTerms/ERF_NumericalDiffusion.H b/Source/SourceTerms/ERF_NumericalDiffusion.H index 14212e417..03666b002 100644 --- a/Source/SourceTerms/ERF_NumericalDiffusion.H +++ b/Source/SourceTerms/ERF_NumericalDiffusion.H @@ -5,23 +5,46 @@ #include #include -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& data, - const amrex::Array4< amrex::Real>& rhs, - const amrex::Array4& mf_x, - const amrex::Array4& 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& 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& prim_data, + const amrex::Array4& cell_data, + const amrex::Array4< amrex::Real>& rhs, + const amrex::Array4& mf); + +void NumericalDiffusion_Xmom (const amrex::Box& bx, + const amrex::Real dt, + const amrex::Real num_diff_coeff, + const amrex::Array4& prim_data, + const amrex::Array4& cell_data, + const amrex::Array4< amrex::Real>& rhs, + const amrex::Array4& mf); + +void NumericalDiffusion_Ymom (const amrex::Box& bx, + const amrex::Real dt, + const amrex::Real num_diff_coeff, + const amrex::Array4& prim_data, + const amrex::Array4& cell_data, + const amrex::Array4< amrex::Real>& rhs, + const amrex::Array4& mf); #endif diff --git a/Source/SourceTerms/ERF_NumericalDiffusion.cpp b/Source/SourceTerms/ERF_NumericalDiffusion.cpp index e5d287d40..5e1bb4583 100644 --- a/Source/SourceTerms/ERF_NumericalDiffusion.cpp +++ b/Source/SourceTerms/ERF_NumericalDiffusion.cpp @@ -9,76 +9,130 @@ 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& data, - const Array4< Real>& rhs, - const Array4& mf_x, - const Array4& 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& prim_data, + const Array4& cell_data, + const Array4< Real>& rhs, + const Array4& 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& mfx_arr = mf_x_bar.array(); - const Array4& 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 * ( (xflux_hi - xflux_lo) * mf_arr(i,j,0) + + (yflux_hi - yflux_lo) * mf_arr(i,j,0) ); }); +} + +/** + * 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& prim_data, + const Array4& cell_data, + const Array4< Real>& rhs, + const Array4& 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 * ( (xflux_hi - xflux_lo) * mf_arr(i,j,0) + + (yflux_hi - yflux_lo) * mf_arr(i,j,0) ); }); } + /** * Function to compute 6th order numerical diffusion RHS. * @@ -86,35 +140,59 @@ NumericalDiffusion (const Box& bx, * @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& data, - const Array4< Real>& rhs) +NumericalDiffusion_Ymom (const Box& bx, + const Real dt, + const Real num_diff_coeff, + const Array4& prim_data, + const Array4& cell_data, + const Array4< Real>& rhs, + const Array4& 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 * ( (xflux_hi - xflux_lo) * mf_arr(i,j,0) + + (yflux_hi - yflux_lo) * mf_arr(i,j,0) ); }); } diff --git a/Source/SourceTerms/ERF_Src_headers.H b/Source/SourceTerms/ERF_Src_headers.H index 7d19eeb0f..667a9dfe7 100644 --- a/Source/SourceTerms/ERF_Src_headers.H +++ b/Source/SourceTerms/ERF_Src_headers.H @@ -37,6 +37,7 @@ void make_sources (int level, int nrk, const SolverChoice& solverChoice, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, + std::unique_ptr& mapfac_m, const amrex::Real* dptr_rhotheta_src, const amrex::Real* dptr_rhoqt_src, const amrex::Real* dptr_wbar_sub, diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index 2e9a4d31f..76c80fb81 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -43,15 +43,15 @@ void make_mom_sources (int level, const MultiFab & S_prim, std::unique_ptr& z_phys_nd, std::unique_ptr& z_phys_cc, - const MultiFab & /*xvel*/, - const MultiFab & /*yvel*/, + const MultiFab & xvel, + const MultiFab & yvel, MultiFab & xmom_src, MultiFab & ymom_src, MultiFab & zmom_src, const MultiFab& base_state, const Geometry geom, const SolverChoice& solverChoice, - std::unique_ptr& mapfac_m, + std::unique_ptr& /*mapfac_m*/, std::unique_ptr& mapfac_u, std::unique_ptr& mapfac_v, const Real* dptr_u_geos, @@ -208,12 +208,15 @@ void make_mom_sources (int level, const Array4& rho_v = S_data[IntVars::ymom].array(mfi); const Array4& rho_w = S_data[IntVars::zmom].array(mfi); + const Array4& u = xvel.array(mfi); + const Array4& v = yvel.array(mfi); + const Array4< Real>& xmom_src_arr = xmom_src.array(mfi); const Array4< Real>& ymom_src_arr = ymom_src.array(mfi); const Array4< Real>& zmom_src_arr = zmom_src.array(mfi); // Map factors - const Array4& mf_m = mapfac_m->const_array(mfi); + //const Array4& mf_m = mapfac_m->const_array(mfi); const Array4& mf_u = mapfac_u->const_array(mfi); const Array4& mf_v = mapfac_v->const_array(mfi); @@ -424,12 +427,12 @@ void make_mom_sources (int level, // 7. Add NUMERICAL DIFFUSION terms // ***************************************************************************** if (l_use_ndiff) { - NumericalDiffusion(tbx, 0, 1, dt, solverChoice.NumDiffCoeff, - rho_u, xmom_src_arr, mf_m, mf_v, false, true); - NumericalDiffusion(tby, 0, 1, dt, solverChoice.NumDiffCoeff, - rho_v, ymom_src_arr, mf_u, mf_m, true, false); - NumericalDiffusion(tbz, 0, 1, dt, solverChoice.NumDiffCoeff, - rho_w, zmom_src_arr, mf_u, mf_v, false, false); + /* + NumericalDiffusion_Xmom(tbx, dt, solverChoice.NumDiffCoeff, + u, cell_data, xmom_src_arr, mf_u); + NumericalDiffusion_Ymom(tby, dt, solverChoice.NumDiffCoeff, + v, cell_data, ymom_src_arr, mf_v); + */ } // ***************************************************************************** diff --git a/Source/SourceTerms/ERF_make_sources.cpp b/Source/SourceTerms/ERF_make_sources.cpp index 2612849c2..d97cd41c1 100644 --- a/Source/SourceTerms/ERF_make_sources.cpp +++ b/Source/SourceTerms/ERF_make_sources.cpp @@ -40,8 +40,9 @@ void make_sources (int level, #endif const Geometry geom, const SolverChoice& solverChoice, - std::unique_ptr& mapfac_u, - std::unique_ptr& mapfac_v, + std::unique_ptr& /*mapfac_u*/, + std::unique_ptr& /*mapfac_v*/, + std::unique_ptr& mapfac_m, const Real* dptr_rhotheta_src, const Real* dptr_rhoqt_src, const Real* dptr_wbar_sub, @@ -330,27 +331,27 @@ void make_sources (int level, // 6. Add numerical diffuion for rho and (rho theta) // ************************************************************************************* if (l_use_ndiff) { - int start_comp = 0; - int num_comp = 2; + int sc; + int nc; - const Array4& mf_u = mapfac_u->const_array(mfi); - const Array4& mf_v = mapfac_v->const_array(mfi); + const Array4& mf_m = mapfac_m->const_array(mfi); + + // Rho is a special case + NumericalDiffusion_Scal(bx, sc=0, nc=1, dt, solverChoice.NumDiffCoeff, + cell_data, cell_data, cell_src, mf_m); + + // Other scalars proceed as normal + NumericalDiffusion_Scal(bx, sc=1, nc=1, dt, solverChoice.NumDiffCoeff, + cell_prim, cell_data, cell_src, mf_m); - NumericalDiffusion(bx, start_comp, num_comp, dt, solverChoice.NumDiffCoeff, - cell_data, cell_src, mf_u, mf_v, false, false); if (l_use_KE && l_diff_KE) { - int sc = RhoKE_comp; - int nc = 1; - NumericalDiffusion(bx, sc, nc, dt, solverChoice.NumDiffCoeff, - cell_data, cell_src, mf_u, mf_v, false, false); - } - { - int sc = RhoScalar_comp; - int nc = NSCALARS; - NumericalDiffusion(bx, sc, nc, dt, solverChoice.NumDiffCoeff, - cell_data, cell_src, mf_u, mf_v, false, false); + NumericalDiffusion_Scal(bx, sc=RhoKE_comp, nc=1, dt, solverChoice.NumDiffCoeff, + cell_prim, cell_data, cell_src, mf_m); } + + NumericalDiffusion_Scal(bx, sc=RhoScalar_comp, nc=NSCALARS, dt, solverChoice.NumDiffCoeff, + cell_prim, cell_data, cell_src, mf_m); } // ************************************************************************************* diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index af489069c..d26784ad7 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -53,7 +53,7 @@ qheating_rates[level].get(), #endif fine_geom, solverChoice, - mapfac_u[level], mapfac_v[level], + mapfac_u[level], mapfac_v[level], mapfac_m[level], dptr_rhotheta_src, dptr_rhoqt_src, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, input_sounding_data, turbPert); @@ -433,7 +433,7 @@ qheating_rates[level], #endif fine_geom, solverChoice, - mapfac_u[level], mapfac_v[level], + mapfac_u[level], mapfac_v[level], mapfac_m[level], dptr_rhotheta_src, dptr_rhoqt_src, dptr_wbar_sub, d_rayleigh_ptrs_at_lev, input_sounding_data, turbPert);