diff --git a/Source/Advection/Advection.H b/Source/Advection/Advection.H index 448c392e8..087ce4ac6 100644 --- a/Source/Advection/Advection.H +++ b/Source/Advection/Advection.H @@ -55,6 +55,7 @@ void AdvectionSrcForMom (const amrex::Box& bxx, const amrex::Box& bxy, const amr const amrex::Array4& w , const amrex::Array4& rho_u , const amrex::Array4& rho_v, const amrex::Array4& Omega , + const amrex::Array4& z_nd, const amrex::Array4& ax, const amrex::Array4& ay, const amrex::Array4& az, diff --git a/Source/Advection/AdvectionSrcForMom.cpp b/Source/Advection/AdvectionSrcForMom.cpp index 0735f66b3..c133405b5 100644 --- a/Source/Advection/AdvectionSrcForMom.cpp +++ b/Source/Advection/AdvectionSrcForMom.cpp @@ -24,6 +24,7 @@ using namespace amrex; * @param[in] rho_u x-component of the momentum * @param[in] rho_v y-component of the momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface + * @param[in] z_nd height coordinate at nodes * @param[in] ax Area fraction of x-faces * @param[in] ay Area fraction of y-faces * @param[in] az Area fraction of z-faces @@ -48,6 +49,7 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, const Array4& rho_u, const Array4& rho_v, const Array4& Omega, + const Array4& z_nd, const Array4& ax, const Array4& ay, const Array4& az, @@ -215,11 +217,13 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, Real xflux_lo = 0.25 * (rho_u(i,j,k) * mf_u_inv(i,j,0) + rho_u(i-1,j,k) * mf_u_inv(i-1,j,0)) * (u(i-1,j,k) + u(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i-1,j,k)); - Real yflux_hi = 0.25 * (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)) * - (u(i,j+1,k) + u(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j+1,k)); + Real met_h_zeta_yhi = Compute_h_zeta_AtEdgeCenterK(i,j+1,k,cellSizeInv,z_nd); + Real yflux_hi = 0.25 * (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)) * + (u(i,j+1,k) + u(i,j,k)) * met_h_zeta_yhi; - Real yflux_lo = 0.25 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i-1,j,k) * mf_v_inv(i-1,j,0)) * - (u(i,j-1,k) + u(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j-1,k)); + Real met_h_zeta_ylo = Compute_h_zeta_AtEdgeCenterK(i,j ,k,cellSizeInv,z_nd); + Real yflux_lo = 0.25 * (rho_v(i,j ,k)*mf_v_inv(i,j ,0) + rho_v(i-1,j ,k)*mf_v_inv(i-1,j ,0)) * + (u(i,j-1,k) + u(i,j,k)) * met_h_zeta_ylo; Real zflux_hi = 0.25 * (Omega(i,j,k+1) + Omega(i-1,j,k+1)) * (u(i,j,k+1) + u(i,j,k)) * 0.5 * (az(i,j,k+1) + az(i-1,j,k+1)); @@ -236,16 +240,19 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xflux_hi = 0.25 * (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)) * - (v(i+1,j,k) + v(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i+1,j,k)); - Real xflux_lo = 0.25 * (rho_u(i,j,k) * mf_u_inv(i,j,0) + rho_u(i ,j-1, k) * mf_u_inv(i-1,j ,0)) * - (v(i-1,j,k) + v(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i-1,j,k)); + Real met_h_zeta_xhi = Compute_h_zeta_AtEdgeCenterK(i+1,j,k,cellSizeInv,z_nd); + Real xflux_hi = 0.25 * (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)) * + (v(i+1,j,k) + v(i,j,k)) * met_h_zeta_xhi; - Real yflux_hi = 0.25 * (rho_v(i,j+1,k) * mf_v_inv(i,j+1,0) + rho_v(i ,j ,k) * mf_v_inv(i ,j ,0)) * + Real met_h_zeta_xlo = Compute_h_zeta_AtEdgeCenterK(i ,j,k,cellSizeInv,z_nd); + Real xflux_lo = 0.25 * (rho_u(i, j, k)*mf_u_inv(i ,j,0) + rho_u(i ,j-1,k)*mf_u_inv(i-1,j ,0)) * + (v(i-1,j,k) + v(i,j,k)) * met_h_zeta_xlo; + + Real yflux_hi = 0.25 * (rho_v(i,j+1,k)*mf_v_inv(i,j+1,0) + rho_v(i,j ,k) * mf_v_inv(i,j ,0)) * (v(i,j+1,k) + v(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j+1,k)); - Real yflux_lo = 0.25 * (rho_v(i,j,k) * mf_v_inv(i,j,0) + rho_v(i , j-1, k) * mf_v_inv(i ,j-1,0)) * + Real yflux_lo = 0.25 * (rho_v(i,j ,k)*mf_v_inv(i,j ,0) + rho_v(i,j-1,k) * mf_v_inv(i,j-1,0)) * (v(i,j-1,k) + v(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j-1,k)); Real zflux_hi = 0.25 * (Omega(i,j,k+1) + Omega(i, j-1, k+1)) * (v(i,j,k+1) + v(i,j,k)) * @@ -263,17 +270,21 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - Real xflux_hi = 0.25*(rho_u(i+1,j,k) + rho_u(i+1,j,k-1)) * mf_u_inv(i+1,j ,0) * - (w(i+1,j,k) + w(i,j,k)) * 0.5 * (ax(i+1,j,k) + ax(i+1,j,k-1)); + Real met_h_zeta_xhi = Compute_h_zeta_AtEdgeCenterJ(i+1,j ,k ,cellSizeInv,z_nd); + Real xflux_hi = 0.25*(rho_u(i+1,j ,k) + rho_u(i+1,j,k-1)) * mf_u_inv(i+1,j,0) * + (w(i+1,j,k) + w(i,j,k)) * met_h_zeta_xhi; - Real xflux_lo = 0.25*(rho_u(i,j,k) + rho_u(i,j,k-1)) * mf_u_inv(i ,j ,0) * - (w(i-1,j,k) + w(i,j,k)) * 0.5 * (ax(i,j,k) + ax(i,j,k-1)); + Real met_h_zeta_xlo = Compute_h_zeta_AtEdgeCenterJ(i ,j ,k ,cellSizeInv,z_nd); + Real xflux_lo = 0.25*(rho_u(i ,j ,k) + rho_u(i ,j,k-1)) * mf_u_inv(i ,j,0) * + (w(i-1,j,k) + w(i,j,k)) * met_h_zeta_xlo; - Real yflux_hi = 0.25*(rho_v(i ,j+1,k) + rho_v(i, j+1, k-1)) * mf_v_inv(i ,j+1,0) * - (w(i,j+1,k) + w(i,j,k)) * 0.5 * (ay(i,j+1,k) + ay(i,j+1,k-1)); + Real met_h_zeta_yhi = Compute_h_zeta_AtEdgeCenterI(i ,j+1,k ,cellSizeInv,z_nd); + Real yflux_hi = 0.25*(rho_v(i,j+1,k) + rho_v(i,j+1,k-1)) * mf_v_inv(i,j+1,0) * + (w(i,j+1,k) + w(i,j,k)) * met_h_zeta_yhi; - Real yflux_lo = 0.25*(rho_v(i ,j ,k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0) * - (w(i,j-1,k) + w(i,j,k)) * 0.5 * (ay(i,j,k) + ay(i,j,k-1)); + Real met_h_zeta_ylo = Compute_h_zeta_AtEdgeCenterI(i ,j ,k ,cellSizeInv,z_nd); + Real yflux_lo = 0.25*(rho_v(i,j ,k) + rho_v(i,j ,k-1)) * mf_v_inv(i,j ,0) * + (w(i,j-1,k) + w(i,j,k)) * met_h_zeta_ylo; Real zflux_lo = 0.25 * (Omega(i,j,k) + Omega(i,j,k-1)) * (w(i,j,k) + w(i,j,k-1)); @@ -294,35 +305,35 @@ AdvectionSrcForMom (const Box& bxx, const Box& bxy, const Box& bxz, if (horiz_adv_type == AdvType::Centered_2nd) { AdvectionSrcForMomVert(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, 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(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, 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::Centered_4th) { AdvectionSrcForMomVert(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, 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(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, 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::Centered_6th) { AdvectionSrcForMomVert(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, horiz_upw_frac, vert_upw_frac, vert_adv_type, lo_z_face, hi_z_face); diff --git a/Source/Advection/AdvectionSrcForMom_N.H b/Source/Advection/AdvectionSrcForMom_N.H index c2fbec130..b359f6777 100644 --- a/Source/Advection/AdvectionSrcForMom_N.H +++ b/Source/Advection/AdvectionSrcForMom_N.H @@ -148,7 +148,7 @@ AdvectionSrcForYMom_N (int i, int j, int k, + (yflux_hi - yflux_lo) * dyInv * mfsq + (zflux_hi - zflux_lo) * dzInv; - return advectionSrc; + return advectionSrc; } /** @@ -270,7 +270,7 @@ AdvectionSrcForZMom_N (int i, int j, int k, + (yflux_hi - yflux_lo) * dyInv * mfsq + (zflux_hi - zflux_lo) * dzInv; - return advectionSrc; + return advectionSrc; } /** diff --git a/Source/Advection/AdvectionSrcForMom_T.H b/Source/Advection/AdvectionSrcForMom_T.H index fdbbf6fca..a16a2129c 100644 --- a/Source/Advection/AdvectionSrcForMom_T.H +++ b/Source/Advection/AdvectionSrcForMom_T.H @@ -10,7 +10,7 @@ * @param[in] rho_u x-component of momentum * @param[in] rho_v y-component of momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface - * @param[in] u x-component of velocity + * @param[in] z_nd height coordinate at nodes * @param[in] ax Area fractions on x-faces * @param[in] ay Area fractions on y-faces * @param[in] az Area fractions on z-faces @@ -27,6 +27,7 @@ AdvectionSrcForXMom (int i, int j, int k, const amrex::Array4& rho_u, const amrex::Array4& rho_v, const amrex::Array4& Omega, + const amrex::Array4& z_nd, const amrex::Array4& ax, const amrex::Array4& ay, const amrex::Array4& az, @@ -70,8 +71,8 @@ AdvectionSrcForXMom (int i, int j, int k, 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); - amrex::Real edgeFluxXYNext = rho_v_avg_hi * interp_hi * 0.5 * (ay(i,j+1,k) + ay(i-1,j+1,k)); - amrex::Real edgeFluxXYPrev = rho_v_avg_lo * interp_lo * 0.5 * (ay(i,j ,k) + ay(i-1,j ,k)); + amrex::Real edgeFluxXYNext = rho_v_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterK(i,j+1,k,cellSizeInv,z_nd); + amrex::Real edgeFluxXYPrev = rho_v_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterK(i,j ,k,cellSizeInv,z_nd); // **************************************************************************************** // Z-fluxes (at edges in j-direction) @@ -105,7 +106,7 @@ AdvectionSrcForXMom (int i, int j, int k, * @param[in] rho_u x-component of momentum * @param[in] rho_v y-component of momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface - * @param[in] v y-component of velocity + * @param[in] z_nd height coordinate at nodes * @param[in] ax Area fractions on x-faces * @param[in] ay Area fractions on y-faces * @param[in] az Area fractions on z-faces @@ -122,6 +123,7 @@ AdvectionSrcForYMom (int i, int j, int k, const amrex::Array4& rho_u, const amrex::Array4& rho_v, const amrex::Array4& Omega, + const amrex::Array4& z_nd, const amrex::Array4& ax, const amrex::Array4& ay, const amrex::Array4& az, @@ -152,8 +154,8 @@ AdvectionSrcForYMom (int i, int j, int k, interp_v_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h); interp_v_h.InterpolateInX(i ,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h); - amrex::Real edgeFluxYXNext = rho_u_avg_hi * 0.5 * (ax(i+1,j,k) + ax(i+1,j-1,k)) * interp_hi; - amrex::Real edgeFluxYXPrev = rho_u_avg_lo * 0.5 * (ax(i ,j,k) + ax(i ,j-1,k)) * interp_lo; + amrex::Real edgeFluxYXNext = rho_u_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterK(i+1,j,k,cellSizeInv,z_nd); + amrex::Real edgeFluxYXPrev = rho_u_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterK(i ,j,k,cellSizeInv,z_nd); // **************************************************************************************** // y-fluxes (at cell centers) @@ -200,6 +202,7 @@ AdvectionSrcForYMom (int i, int j, int k, * @param[in] rho_v y-component of momentum * @param[in] Omega component of the momentum normal to the z-coordinate surface * @param[in] w z-component of velocity + * @param[in] z_nd height coordinate at nodes * @param[in] ax Area fractions on x-faces * @param[in] ay Area fractions on y-faces * @param[in] az Area fractions on z-faces @@ -221,6 +224,7 @@ AdvectionSrcForZMom (int i, int j, int k, const amrex::Array4& rho_v, const amrex::Array4& Omega, const amrex::Array4& w, + const amrex::Array4& z_nd, const amrex::Array4& ax, const amrex::Array4& ay, const amrex::Array4& az, @@ -255,8 +259,8 @@ AdvectionSrcForZMom (int i, int j, int k, interp_omega_h.InterpolateInX(i+1,j,k,0,interp_hi,rho_u_avg_hi,upw_frac_h); interp_omega_h.InterpolateInX(i,j,k,0,interp_lo,rho_u_avg_lo,upw_frac_h); - amrex::Real edgeFluxZXNext = rho_u_avg_hi * 0.5 * (ax(i+1,j,k) + ax(i+1,j,k-1)) * interp_hi; - amrex::Real edgeFluxZXPrev = rho_u_avg_lo * 0.5 * (ax(i ,j,k) + ax(i ,j,k-1)) * interp_lo; + amrex::Real edgeFluxZXNext = rho_u_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterJ(i+1,j,k,cellSizeInv,z_nd); + amrex::Real edgeFluxZXPrev = rho_u_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterJ(i ,j,k,cellSizeInv,z_nd); // **************************************************************************************** // y-fluxes (at edges in i-direction) @@ -267,8 +271,8 @@ AdvectionSrcForZMom (int i, int j, int k, rho_v_avg_lo = 0.5 * (rho_v(i, j , k) + rho_v(i, j , k-1)) * mf_v_inv(i ,j ,0); interp_omega_h.InterpolateInY(i,j,k,0,interp_lo,rho_v_avg_lo,upw_frac_h); - amrex::Real edgeFluxZYNext = rho_v_avg_hi * 0.5 * (ay(i,j+1,k) + ay(i,j+1,k-1)) * interp_hi; - amrex::Real edgeFluxZYPrev = rho_v_avg_lo * 0.5 * (ay(i,j ,k) + ay(i,j ,k-1)) * interp_lo; + amrex::Real edgeFluxZYNext = rho_v_avg_hi * interp_hi * Compute_h_zeta_AtEdgeCenterI(i,j+1,k,cellSizeInv,z_nd); + amrex::Real edgeFluxZYPrev = rho_v_avg_lo * interp_lo * Compute_h_zeta_AtEdgeCenterI(i,j ,k,cellSizeInv,z_nd); // **************************************************************************************** // z-fluxes (at cell centers) @@ -356,6 +360,7 @@ AdvectionSrcForMomWrapper (const amrex::Box& bxx, const amrex::Box& bxy, const a const amrex::Array4& u, const amrex::Array4& v, const amrex::Array4& w, + const amrex::Array4& z_nd, const amrex::Array4& ax, const amrex::Array4& ay, const amrex::Array4& az, @@ -378,21 +383,21 @@ AdvectionSrcForMomWrapper (const amrex::Box& bxx, const amrex::Box& bxy, const a amrex::ParallelFor(bxx, bxy, bxz, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - rho_u_rhs(i, j, k) = -AdvectionSrcForXMom(i, j, k, rho_u, rho_v, Omega, ax, ay, az, detJ, + rho_u_rhs(i, j, k) = -AdvectionSrcForXMom(i, j, k, rho_u, rho_v, Omega, z_nd, ax, ay, az, detJ, interp_u_h, interp_u_v, upw_frac_h, upw_frac_v, cellSizeInv, mf_u_inv, mf_v_inv); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - rho_v_rhs(i, j, k) = -AdvectionSrcForYMom(i, j, k, rho_u, rho_v, Omega, ax, ay, az, detJ, + rho_v_rhs(i, j, k) = -AdvectionSrcForYMom(i, j, k, rho_u, rho_v, Omega, z_nd, ax, ay, az, detJ, interp_v_h, interp_v_v, upw_frac_h, upw_frac_v, cellSizeInv, mf_u_inv, mf_v_inv); }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - rho_w_rhs(i, j, k) = -AdvectionSrcForZMom(i, j, k, rho_u, rho_v, Omega, w, ax, ay, az, detJ, + rho_w_rhs(i, j, k) = -AdvectionSrcForZMom(i, j, k, rho_u, rho_v, Omega, w, z_nd, ax, ay, az, detJ, interp_w_h, interp_w_v, interp_w_wall, upw_frac_h, upw_frac_v, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, @@ -415,6 +420,7 @@ AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amre const amrex::Array4& u, const amrex::Array4& v, const amrex::Array4& w, + const amrex::Array4& z_nd, const amrex::Array4& ax, const amrex::Array4& ay, const amrex::Array4& az, @@ -431,7 +437,7 @@ AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amre if (vert_adv_type == AdvType::Centered_2nd) { AdvectionSrcForMomWrapper(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, upw_frac_h, upw_frac_v, @@ -439,7 +445,7 @@ AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amre } else if (vert_adv_type == AdvType::Upwind_3rd) { AdvectionSrcForMomWrapper(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, upw_frac_h, upw_frac_v, @@ -447,7 +453,7 @@ AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amre } else if (vert_adv_type == AdvType::Centered_4th) { AdvectionSrcForMomWrapper(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, upw_frac_h, upw_frac_v, @@ -455,7 +461,7 @@ AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amre } else if (vert_adv_type == AdvType::Upwind_5th) { AdvectionSrcForMomWrapper(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, upw_frac_h, upw_frac_v, @@ -463,7 +469,7 @@ AdvectionSrcForMomVert (const amrex::Box& bxx, const amrex::Box& bxy, const amre } else if (vert_adv_type == AdvType::Centered_6th) { AdvectionSrcForMomWrapper(bxx, bxy, bxz, rho_u_rhs, rho_v_rhs, rho_w_rhs, - rho_u, rho_v, Omega, u, v, w, ax, ay, az, detJ, + rho_u, rho_v, Omega, u, v, w, z_nd, ax, ay, az, detJ, cellSizeInv, mf_m, mf_u_inv, mf_v_inv, upw_frac_h, upw_frac_v, diff --git a/Source/Diffusion/ComputeStrain_T.cpp b/Source/Diffusion/ComputeStrain_T.cpp index bc069d232..bc22b56a4 100644 --- a/Source/Diffusion/ComputeStrain_T.cpp +++ b/Source/Diffusion/ComputeStrain_T.cpp @@ -234,9 +234,9 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, Real GradWz = 0.25 * dxInv[2] * ( w(i ,j ,k+1) + w(i ,j-1,k+1) - w(i ,j ,k-1) - w(i ,j-1,k-1) ); - Real met_h_eta,met_h_zeta; - met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd); - met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd); + Real met_h_eta,met_h_zeta; + met_h_eta = Compute_h_eta_AtEdgeCenterI (i,j,k,dxInv,z_nd); + met_h_zeta = Compute_h_zeta_AtEdgeCenterI(i,j,k,dxInv,z_nd); tau23(i,j,k) = 0.5 * ( (v(i, j, k) - v(i, j, k-1))*dxInv[2]/met_h_zeta + ( (-(8./3.) * w(i,j-1,k) + 3. * w(i,j ,k) - (1./3.) * w(i,j+1,k))*dxInv[1]*mf_v(i,j,0) @@ -358,8 +358,7 @@ ComputeStrain_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Box domain, Real met_h_xi,met_h_eta,met_h_zeta; met_h_xi = Compute_h_xi_AtEdgeCenterK (i,j,k,dxInv,z_nd); met_h_eta = Compute_h_eta_AtEdgeCenterK (i,j,k,dxInv,z_nd); - // met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); - met_h_zeta = 0.5 * (detJ(i,j,k) + detJ(i,j,k-1)); + met_h_zeta = Compute_h_zeta_AtEdgeCenterK(i,j,k,dxInv,z_nd); tau12(i,j,k) = 0.5 * ( (u(i, j, k)/mf_u(i,j,0) - u(i , j-1, k)/mf_u(i,j-1,0))*dxInv[1] + (v(i, j, k)/mf_v(i,j,0) - v(i-1, j , k)/mf_v(i-1,j,0))*dxInv[0] diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 9db14dc2d..cb9daf933 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -838,7 +838,7 @@ void erf_slow_rhs_pre (int level, int finest_level, rho_u_rhs, rho_v_rhs, rho_w_rhs, cell_data, u, v, w, rho_u, rho_v, omega_arr, - ax_arr, ay_arr, az_arr, detJ_arr, + z_nd, ax_arr, ay_arr, az_arr, detJ_arr, dxInv, mf_m, mf_u, mf_v, l_horiz_adv_type, l_vert_adv_type, l_horiz_upw_frac, l_vert_upw_frac,