diff --git a/Source/Diffusion/ERF_ComputeStress_N.cpp b/Source/Diffusion/ERF_ComputeStress_N.cpp index 1eef5b552..7985a98e1 100644 --- a/Source/Diffusion/ERF_ComputeStress_N.cpp +++ b/Source/Diffusion/ERF_ComputeStress_N.cpp @@ -100,6 +100,7 @@ ComputeStressConsVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, * @param[in,out] tau13 13 strain -> stress * @param[in,out] tau23 23 strain -> stress * @param[in] er_arr expansion rate + * @param[in] horizontal_only ignore vertical viscous stresses */ void ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, @@ -107,7 +108,7 @@ ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, const Array4& cell_data, Array4& tau11, Array4& tau22, Array4& tau33, Array4& tau12, Array4& tau13, Array4& tau23, - const Array4& er_arr) + const Array4& er_arr, const bool horizontal_only) { Real OneThird = (1./3.); @@ -119,7 +120,7 @@ ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, Real rhoAlpha = cell_data(i, j, k, Rho_comp) * mu_eff; Real mu_11 = rhoAlpha + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_h); Real mu_22 = mu_11; - Real mu_33 = rhoAlpha + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_v); + Real mu_33 = horizontal_only ? 0.0 : rhoAlpha + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_v); tau11(i,j,k) = -mu_11 * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) ); tau22(i,j,k) = -mu_22 * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) ); tau33(i,j,k) = -mu_33 * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) ); @@ -140,7 +141,7 @@ ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, + cell_data(i-1, j, k-1, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real mu_13 = rho_bar*mu_eff + 2.0*mu_bar; + Real mu_13 = horizontal_only ? 0.0 : rho_bar*mu_eff + 2.0*mu_bar; tau13(i,j,k) *= -mu_13; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { @@ -148,7 +149,7 @@ ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, + cell_data(i, j-1, k-1, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real mu_23 = rho_bar*mu_eff + 2.0*mu_bar; + Real mu_23 = horizontal_only ? 0.0 : rho_bar*mu_eff + 2.0*mu_bar; tau23(i,j,k) *= -mu_23; }); } @@ -159,7 +160,7 @@ ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, ParallelFor(bxcc, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real mu_11 = mu_eff + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_h); Real mu_22 = mu_11; - Real mu_33 = mu_eff + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_v); + Real mu_33 = horizontal_only ? 0.0 : mu_eff + 2.0 * mu_turb(i, j, k, EddyDiff::Mom_v); tau11(i,j,k) = -mu_11 * ( tau11(i,j,k) - OneThird*er_arr(i,j,k) ); tau22(i,j,k) = -mu_22 * ( tau22(i,j,k) - OneThird*er_arr(i,j,k) ); tau33(i,j,k) = -mu_33 * ( tau33(i,j,k) - OneThird*er_arr(i,j,k) ); @@ -176,13 +177,13 @@ ComputeStressVarVisc_N (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real mu_13 = mu_eff + 2.0*mu_bar; + Real mu_13 = horizontal_only ? 0.0 : mu_eff + 2.0*mu_bar; tau13(i,j,k) *= -mu_13; }, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real mu_23 = mu_eff + 2.0*mu_bar; + Real mu_23 = horizontal_only ? 0.0 : mu_eff + 2.0*mu_bar; tau23(i,j,k) *= -mu_23; }); } diff --git a/Source/Diffusion/ERF_ComputeStress_T.cpp b/Source/Diffusion/ERF_ComputeStress_T.cpp index 391cc47f8..02b8eb3a4 100644 --- a/Source/Diffusion/ERF_ComputeStress_T.cpp +++ b/Source/Diffusion/ERF_ComputeStress_T.cpp @@ -298,6 +298,7 @@ ComputeStressConsVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, * @param[in] er_arr expansion rate * @param[in] z_nd nodal array of physical z heights * @param[in] dxInv inverse cell size array + * @param[in] horizontal_only ignore vertical viscous stresses */ void ComputeStressVarVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, @@ -310,7 +311,7 @@ ComputeStressVarVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, const Array4& er_arr, const Array4& z_nd, const Array4& detJ, - const GpuArray& dxInv) + const GpuArray& dxInv, const bool horizontal_only) { // Handle constant alpha case, in which the provided mu_eff is actually // "alpha" and the viscosity needs to be scaled by rho. This can be further @@ -350,205 +351,223 @@ ComputeStressVarVisc_T (Box bxcc, Box tbxxy, Box tbxxz, Box tbxyz, Real mu_eff, tau33(i,j,k) -= OneThird*er_arr(i,j,k); }); - // Second block: compute 2mu*JT*(S-D) - //*********************************************************************************** - // Fill tau33 first (no linear combination extrapolation) - //----------------------------------------------------------------------------------- - ParallelFor(bxcc2, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_xi,met_h_eta; - met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd); - met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd); - - Real tau31bar = 0.25 * ( tau31(i , j , k ) + tau31(i+1, j , k ) - + tau31(i , j , k+1) + tau31(i+1, j , k+1) ); - Real tau32bar = 0.25 * ( tau32(i , j , k ) + tau32(i , j+1, k ) - + tau32(i , j , k+1) + tau32(i , j+1, k+1) ); - - Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_v); - - tau33(i,j,k) -= met_h_xi*tau31bar + met_h_eta*tau32bar; - tau33(i,j,k) *= -mu_tot; - }); - - // Second block: compute 2mu*JT*(S-D) - //*********************************************************************************** - // Fill tau13, tau23 next (linear combination extrapolation) - //----------------------------------------------------------------------------------- - // Extrapolate tau13 & tau23 to bottom - { - Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) ); - tbxxz.growLo(2,-1); - ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + if (!horizontal_only) { + // Second block: compute 2mu*JT*(S-D) + //*********************************************************************************** + // Fill tau33 first (no linear combination extrapolation) + //----------------------------------------------------------------------------------- + ParallelFor(bxcc2, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real met_h_xi,met_h_eta; + met_h_xi = Compute_h_xi_AtCellCenter (i,j,k,dxInv,z_nd); + met_h_eta = Compute_h_eta_AtCellCenter (i,j,k,dxInv,z_nd); + + Real tau31bar = 0.25 * ( tau31(i , j , k ) + tau31(i+1, j , k ) + + tau31(i , j , k+1) + tau31(i+1, j , k+1) ); + Real tau32bar = 0.25 * ( tau32(i , j , k ) + tau32(i , j+1, k ) + + tau32(i , j , k+1) + tau32(i , j+1, k+1) ); + + Real mu_tot = rhoAlpha(i,j,k) + 2.0*mu_turb(i, j, k, EddyDiff::Mom_v); + + tau33(i,j,k) -= met_h_xi*tau31bar + met_h_eta*tau32bar; + tau33(i,j,k) *= -mu_tot; + }); + + // Second block: compute 2mu*JT*(S-D) + //*********************************************************************************** + // Fill tau13, tau23 next (linear combination extrapolation) + //----------------------------------------------------------------------------------- + // Extrapolate tau13 & tau23 to bottom { - Real met_h_xi,met_h_eta,met_h_zeta; - met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd); - met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd); - met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd); - - Real tau11lo = 0.5 * ( tau11(i , j , k ) + tau11(i-1, j , k ) ); - Real tau11hi = 0.5 * ( tau11(i , j , k+1) + tau11(i-1, j , k+1) ); - Real tau11bar = 1.5*tau11lo - 0.5*tau11hi; - - Real tau12lo = 0.5 * ( tau12(i , j , k ) + tau12(i , j+1, k ) ); - Real tau12hi = 0.5 * ( tau12(i , j , k+1) + tau12(i , j+1, k+1) ); - Real tau12bar = 1.5*tau12lo - 0.5*tau12hi; - - Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) - + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k ) - + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) ); - Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; - - tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar; - tau13(i,j,k) *= -mu_tot; - - tau31(i,j,k) *= -mu_tot*met_h_zeta; - }); - - Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) ); - tbxyz.growLo(2,-1); - ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_xi,met_h_eta,met_h_zeta; - met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd); - 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 tau21lo = 0.5 * ( tau21(i , j , k ) + tau21(i+1, j , k ) ); - Real tau21hi = 0.5 * ( tau21(i , j , k+1) + tau21(i+1, j , k+1) ); - Real tau21bar = 1.5*tau21lo - 0.5*tau21hi; - - Real tau22lo = 0.5 * ( tau22(i , j , k ) + tau22(i , j-1, k ) ); - Real tau22hi = 0.5 * ( tau22(i , j , k+1) + tau22(i , j-1, k+1) ); - Real tau22bar = 1.5*tau22lo - 0.5*tau22hi; - - Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) - + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k ) - + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) ); - Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; - - tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar; - tau23(i,j,k) *= -mu_tot; - - tau32(i,j,k) *= -mu_tot*met_h_zeta; - }); - } - // Extrapolate tau13 & tau23 to top - { - Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) ); - tbxxz.growHi(2,-1); - ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_xi,met_h_eta,met_h_zeta; - met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd); - met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd); - met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd); - - Real tau11lo = 0.5 * ( tau11(i , j , k-2) + tau11(i-1, j , k-2) ); - Real tau11hi = 0.5 * ( tau11(i , j , k-1) + tau11(i-1, j , k-1) ); - Real tau11bar = 1.5*tau11hi - 0.5*tau11lo; - - Real tau12lo = 0.5 * ( tau12(i , j , k-2) + tau12(i , j+1, k-2) ); - Real tau12hi = 0.5 * ( tau12(i , j , k-1) + tau12(i , j+1, k-1) ); - Real tau12bar = 1.5*tau12hi - 0.5*tau12lo; - - Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) - + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k ) - + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) ); - Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; - - tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar; - tau13(i,j,k) *= -mu_tot; - - tau31(i,j,k) *= -mu_tot*met_h_zeta; - }); - - Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) ); - tbxyz.growHi(2,-1); - ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + Box planexz = tbxxz; planexz.setBig(2, planexz.smallEnd(2) ); + tbxxz.growLo(2,-1); + ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real met_h_xi,met_h_eta,met_h_zeta; + met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd); + met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd); + met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd); + + Real tau11lo = 0.5 * ( tau11(i , j , k ) + tau11(i-1, j , k ) ); + Real tau11hi = 0.5 * ( tau11(i , j , k+1) + tau11(i-1, j , k+1) ); + Real tau11bar = 1.5*tau11lo - 0.5*tau11hi; + + Real tau12lo = 0.5 * ( tau12(i , j , k ) + tau12(i , j+1, k ) ); + Real tau12hi = 0.5 * ( tau12(i , j , k+1) + tau12(i , j+1, k+1) ); + Real tau12bar = 1.5*tau12lo - 0.5*tau12hi; + + Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); + Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k ) + + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) ); + Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; + + tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar; + tau13(i,j,k) *= -mu_tot; + + tau31(i,j,k) *= -mu_tot*met_h_zeta; + }); + + Box planeyz = tbxyz; planeyz.setBig(2, planeyz.smallEnd(2) ); + tbxyz.growLo(2,-1); + ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real met_h_xi,met_h_eta,met_h_zeta; + met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd); + 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 tau21lo = 0.5 * ( tau21(i , j , k ) + tau21(i+1, j , k ) ); + Real tau21hi = 0.5 * ( tau21(i , j , k+1) + tau21(i+1, j , k+1) ); + Real tau21bar = 1.5*tau21lo - 0.5*tau21hi; + + Real tau22lo = 0.5 * ( tau22(i , j , k ) + tau22(i , j-1, k ) ); + Real tau22hi = 0.5 * ( tau22(i , j , k+1) + tau22(i , j-1, k+1) ); + Real tau22bar = 1.5*tau22lo - 0.5*tau22hi; + + Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); + Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k ) + + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) ); + Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; + + tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar; + tau23(i,j,k) *= -mu_tot; + + tau32(i,j,k) *= -mu_tot*met_h_zeta; + }); + } + // Extrapolate tau13 & tau23 to top { - Real met_h_xi,met_h_eta,met_h_zeta; - met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd); - 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 tau21lo = 0.5 * ( tau21(i , j , k-2) + tau21(i+1, j , k-2) ); - Real tau21hi = 0.5 * ( tau21(i , j , k-1) + tau21(i+1, j , k-1) ); - Real tau21bar = 1.5*tau21hi - 0.5*tau21lo; - - Real tau22lo = 0.5 * ( tau22(i , j , k-2) + tau22(i , j-1, k-2) ); - Real tau22hi = 0.5 * ( tau22(i , j , k-1) + tau22(i , j-1, k-1) ); - Real tau22bar = 1.5*tau22hi - 0.5*tau22lo; - - Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) - + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); - Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k ) - + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) ); - Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; - - tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar; - tau23(i,j,k) *= -mu_tot; - - tau32(i,j,k) *= -mu_tot*met_h_zeta; - }); + Box planexz = tbxxz; planexz.setSmall(2, planexz.bigEnd(2) ); + tbxxz.growHi(2,-1); + ParallelFor(planexz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real met_h_xi,met_h_eta,met_h_zeta; + met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd); + met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd); + met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd); + + Real tau11lo = 0.5 * ( tau11(i , j , k-2) + tau11(i-1, j , k-2) ); + Real tau11hi = 0.5 * ( tau11(i , j , k-1) + tau11(i-1, j , k-1) ); + Real tau11bar = 1.5*tau11hi - 0.5*tau11lo; + + Real tau12lo = 0.5 * ( tau12(i , j , k-2) + tau12(i , j+1, k-2) ); + Real tau12hi = 0.5 * ( tau12(i , j , k-1) + tau12(i , j+1, k-1) ); + Real tau12bar = 1.5*tau12hi - 0.5*tau12lo; + + Real mu_bar = 0.25*( mu_turb(i-1, j, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + + mu_turb(i-1, j, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); + Real rhoAlpha_bar = 0.25*( rhoAlpha(i-1, j, k ) + rhoAlpha(i, j, k ) + + rhoAlpha(i-1, j, k-1) + rhoAlpha(i, j, k-1) ); + Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; + + tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar; + tau13(i,j,k) *= -mu_tot; + + tau31(i,j,k) *= -mu_tot*met_h_zeta; + }); + + Box planeyz = tbxyz; planeyz.setSmall(2, planeyz.bigEnd(2) ); + tbxyz.growHi(2,-1); + ParallelFor(planeyz,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real met_h_xi,met_h_eta,met_h_zeta; + met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd); + 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 tau21lo = 0.5 * ( tau21(i , j , k-2) + tau21(i+1, j , k-2) ); + Real tau21hi = 0.5 * ( tau21(i , j , k-1) + tau21(i+1, j , k-1) ); + Real tau21bar = 1.5*tau21hi - 0.5*tau21lo; + + Real tau22lo = 0.5 * ( tau22(i , j , k-2) + tau22(i , j-1, k-2) ); + Real tau22hi = 0.5 * ( tau22(i , j , k-1) + tau22(i , j-1, k-1) ); + Real tau22bar = 1.5*tau22hi - 0.5*tau22lo; + + Real mu_bar = 0.25*( mu_turb(i, j-1, k , EddyDiff::Mom_v) + mu_turb(i, j, k , EddyDiff::Mom_v) + + mu_turb(i, j-1, k-1, EddyDiff::Mom_v) + mu_turb(i, j, k-1, EddyDiff::Mom_v) ); + Real rhoAlpha_bar = 0.25*( rhoAlpha(i, j-1, k ) + rhoAlpha(i, j, k ) + + rhoAlpha(i, j-1, k-1) + rhoAlpha(i, j, k-1) ); + Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; + + tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar; + tau23(i,j,k) *= -mu_tot; + + tau32(i,j,k) *= -mu_tot*met_h_zeta; + }); + } + + // Second block: compute 2mu*JT*(S-D) + //*********************************************************************************** + // Fill tau13, tau23 next (valid averaging region) + //----------------------------------------------------------------------------------- + ParallelFor(tbxxz,tbxyz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real met_h_xi,met_h_eta,met_h_zeta; + met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd); + met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd); + met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd); + + Real tau11bar = 0.25 * ( tau11(i , j , k ) + tau11(i-1, j , k ) + + tau11(i , j , k-1) + tau11(i-1, j , k-1) ); + Real tau12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k ) + + tau12(i , j , k-1) + tau12(i , j+1, k-1) ); + + Real mu_bar = 0.25 * ( mu_turb(i-1, j , k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v) + + mu_turb(i-1, j , k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) ); + Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i-1, j , k ) + rhoAlpha(i , j , k ) + + rhoAlpha(i-1, j , k-1) + rhoAlpha(i , j , k-1) ); + Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; + + tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar; + tau13(i,j,k) *= -mu_tot; + + tau31(i,j,k) *= -mu_tot*met_h_zeta; + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real met_h_xi,met_h_eta,met_h_zeta; + met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd); + 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 tau21bar = 0.25 * ( tau21(i , j , k ) + tau21(i+1, j , k ) + + tau21(i , j , k-1) + tau21(i+1, j , k-1) ); + Real tau22bar = 0.25 * ( tau22(i , j , k ) + tau22(i , j-1, k ) + + tau22(i , j , k-1) + tau22(i , j-1, k-1) ); + + Real mu_bar = 0.25 * ( mu_turb(i , j-1, k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v) + + mu_turb(i , j-1, k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) ); + Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i , j-1, k ) + rhoAlpha(i , j , k ) + + rhoAlpha(i , j-1, k-1) + rhoAlpha(i , j , k-1) ); + Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; + + tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar; + tau23(i,j,k) *= -mu_tot; + + tau32(i,j,k) *= -mu_tot*met_h_zeta; + }); + } else { // horizontal only - zero out vertical stresses + ParallelFor(bxcc,tbxxz,tbxyz, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tau33(i,j,k) = 0.0; + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tau13(i,j,k) = 0.0; + tau31(i,j,k) = 0.0; + }, + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + tau23(i,j,k) = 0.0; + tau32(i,j,k) = 0.0; + }); } - // Second block: compute 2mu*JT*(S-D) - //*********************************************************************************** - // Fill tau13, tau23 next (valid averaging region) - //----------------------------------------------------------------------------------- - ParallelFor(tbxxz,tbxyz, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_xi,met_h_eta,met_h_zeta; - met_h_xi = Compute_h_xi_AtEdgeCenterJ (i,j,k,dxInv,z_nd); - met_h_eta = Compute_h_eta_AtEdgeCenterJ (i,j,k,dxInv,z_nd); - met_h_zeta = Compute_h_zeta_AtEdgeCenterJ(i,j,k,dxInv,z_nd); - - Real tau11bar = 0.25 * ( tau11(i , j , k ) + tau11(i-1, j , k ) - + tau11(i , j , k-1) + tau11(i-1, j , k-1) ); - Real tau12bar = 0.25 * ( tau12(i , j , k ) + tau12(i , j+1, k ) - + tau12(i , j , k-1) + tau12(i , j+1, k-1) ); - - Real mu_bar = 0.25 * ( mu_turb(i-1, j , k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v) - + mu_turb(i-1, j , k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) ); - Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i-1, j , k ) + rhoAlpha(i , j , k ) - + rhoAlpha(i-1, j , k-1) + rhoAlpha(i , j , k-1) ); - Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; - - tau13(i,j,k) -= met_h_xi*tau11bar + met_h_eta*tau12bar; - tau13(i,j,k) *= -mu_tot; - - tau31(i,j,k) *= -mu_tot*met_h_zeta; - }, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - Real met_h_xi,met_h_eta,met_h_zeta; - met_h_xi = Compute_h_xi_AtEdgeCenterI (i,j,k,dxInv,z_nd); - 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 tau21bar = 0.25 * ( tau21(i , j , k ) + tau21(i+1, j , k ) - + tau21(i , j , k-1) + tau21(i+1, j , k-1) ); - Real tau22bar = 0.25 * ( tau22(i , j , k ) + tau22(i , j-1, k ) - + tau22(i , j , k-1) + tau22(i , j-1, k-1) ); - - Real mu_bar = 0.25 * ( mu_turb(i , j-1, k , EddyDiff::Mom_v) + mu_turb(i , j , k , EddyDiff::Mom_v) - + mu_turb(i , j-1, k-1, EddyDiff::Mom_v) + mu_turb(i , j , k-1, EddyDiff::Mom_v) ); - Real rhoAlpha_bar = 0.25 * ( rhoAlpha(i , j-1, k ) + rhoAlpha(i , j , k ) - + rhoAlpha(i , j-1, k-1) + rhoAlpha(i , j , k-1) ); - Real mu_tot = rhoAlpha_bar + 2.0*mu_bar; - - tau23(i,j,k) -= met_h_xi*tau21bar + met_h_eta*tau22bar; - tau23(i,j,k) *= -mu_tot; - - tau32(i,j,k) *= -mu_tot*met_h_zeta; - }); - // Fill the remaining components: tau11, tau22, tau12/21 //----------------------------------------------------------------------------------- ParallelFor(bxcc,tbxxy, diff --git a/Source/Diffusion/ERF_Diffusion.H b/Source/Diffusion/ERF_Diffusion.H index e316194d0..b02e24562 100644 --- a/Source/Diffusion/ERF_Diffusion.H +++ b/Source/Diffusion/ERF_Diffusion.H @@ -38,8 +38,6 @@ void DiffusionSrcForMom_T (const amrex::Box& bxx, const amrex::Box& bxy, const a const amrex::Array4& mf_u , const amrex::Array4& mf_v ); - - void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, int start_comp, int num_comp, const bool& exp_most, @@ -66,7 +64,8 @@ void DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& tm_arr, const amrex::GpuArray grav_gpu, const amrex::BCRec* bc_ptr, - const bool use_most); + const bool use_most, + const bool horizontal_only); void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, int start_comp, int num_comp, @@ -104,7 +103,8 @@ void DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, const amrex::Array4& tm_arr, const amrex::GpuArray grav_gpu, const amrex::BCRec* bc_ptr, - const bool use_most); + const bool use_most, + const bool horizontal_only); void ComputeStressConsVisc_N (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxxz, amrex::Box tbxyz, amrex::Real mu_eff, const amrex::Array4& cell_data, @@ -130,7 +130,7 @@ void ComputeStressVarVisc_N (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxxz const amrex::Array4& cell_data, amrex::Array4& tau11, amrex::Array4& tau22, amrex::Array4& tau33, amrex::Array4& tau12, amrex::Array4& tau13, amrex::Array4& tau23, - const amrex::Array4& er_arr); + const amrex::Array4& er_arr, const bool horizontal_only); void ComputeStressVarVisc_T (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxxz, amrex::Box tbxyz, amrex::Real mu_eff, const amrex::Array4& mu_turb, @@ -142,7 +142,7 @@ void ComputeStressVarVisc_T (amrex::Box bxcc, amrex::Box tbxxy, amrex::Box tbxxz const amrex::Array4& er_arr, const amrex::Array4& z_nd, const amrex::Array4& detJ, - const amrex::GpuArray& dxInv); + const amrex::GpuArray& dxInv, const bool horizontal_only); diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp index 58fde87ac..034c77a89 100644 --- a/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/ERF_DiffusionSrcForState_N.cpp @@ -35,6 +35,7 @@ using namespace amrex; * @param[in] grav_gpu gravity vector * @param[in] bc_ptr container with boundary conditions * @param[in] use_most whether we have turned on MOST BCs + * @param[in] horizontal_only don't compute vertical diffusion terms */ void DiffusionSrcForState_N (const Box& bx, const Box& domain, @@ -63,7 +64,8 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, const Array4& tm_arr, const GpuArray grav_gpu, const BCRec* bc_ptr, - const bool use_most) + const bool use_most, + const bool horizontal_only) { BL_PROFILE_VAR("DiffusionSrcForState_N()",DiffusionSrcForState_N); @@ -264,59 +266,61 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; - const int prim_scal_index = prim_index; - - Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; - rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index]) - + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) ); - - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && - k == dom_lo.z); - - if (ext_dir_on_zlo) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + // const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; + const int prim_scal_index = prim_index; + + Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index]) + + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) ); + + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == dom_lo.z); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && + k == dom_lo.z); + + if (ext_dir_on_zlo) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; - } else if (ext_dir_on_zhi) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; - } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; - } - - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; + } else if (ext_dir_on_zhi) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; + } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } - }); + }); + } } else if (l_turb) { // with MolecDiffType::Constant or None ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -381,56 +385,58 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - - Real rhoAlpha = d_alpha_eff[prim_index]; - rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) - + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); - - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && - k == dom_lo.z); - - if (ext_dir_on_zlo) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + + Real rhoAlpha = d_alpha_eff[prim_index]; + rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) + + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); + + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == dom_lo.z); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && + k == dom_lo.z); + + if (ext_dir_on_zlo) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; - } else if (ext_dir_on_zhi) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; - } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; - } - - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; + } else if (ext_dir_on_zhi) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; + } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } - }); + }); + } } else if(l_consA) { // without an LES/PBL model ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -493,55 +499,57 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - - Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && - k == dom_lo.z); - - if (ext_dir_on_zlo) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + + Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; + + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == dom_lo.z); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && + k == dom_lo.z); + + if (ext_dir_on_zlo) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; - } else if (ext_dir_on_zhi) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; - } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; - } - - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; + } else if (ext_dir_on_zhi) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; + } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } - }); + }); + } } else { // with MolecDiffType::Constant or None // without an LES/PBL model @@ -603,53 +611,64 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j-1, k, prim_index)) * dy_inv * mf_v(i,j,0); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - - Real rhoAlpha = d_alpha_eff[prim_index]; - - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == dom_lo.z); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && - k == dom_lo.z); - - if (ext_dir_on_zlo) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + + Real rhoAlpha = d_alpha_eff[prim_index]; + + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == dom_lo.z); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].hi(2) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[BCVars::cons_bc+qty_index].lo(2) == ERFBCType::foextrap) && + k == dom_lo.z); + + if (ext_dir_on_zlo) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; - } else if (ext_dir_on_zhi) { - zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; - } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; - } - - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + - (1./3.) * cell_prim(i, j, k+1, prim_index) ) * dz_inv; + } else if (ext_dir_on_zhi) { + zflux(i,j,k,qty_index) = -rhoAlpha * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ) * dz_inv; + } else if (most_on_zlo && (qty_index == RhoTheta_comp)) { + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * (cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index)) * dz_inv; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } + }); + } + } + + // zero vertical fluxes if doing horizontal only - note we'll deal with MOST flux later for explicit most in this case + if (horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + zflux(i,j,k,qty_index) = 0.0; }); } @@ -659,7 +678,6 @@ DiffusionSrcForState_N (const Box& bx, const Box& domain, int qty_index = start_comp + n; ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { - cell_rhs(i,j,k,qty_index) -= (xflux(i+1,j ,k ,qty_index) - xflux(i, j, k, qty_index)) * dx_inv * mf_m(i,j,0) // Diffusive flux in x-dir +(yflux(i ,j+1,k ,qty_index) - yflux(i, j, k, qty_index)) * dy_inv * mf_m(i,j,0) // Diffusive flux in y-dir +(zflux(i ,j ,k+1,qty_index) - zflux(i, j, k, qty_index)) * dz_inv; // Diffusive flux in z-dir diff --git a/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp b/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp index d210b512c..470a77939 100644 --- a/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/ERF_DiffusionSrcForState_T.cpp @@ -38,6 +38,7 @@ using namespace amrex; * @param[in] grav_gpu gravity vector * @param[in] bc_ptr container with boundary conditions * @param[in] use_most whether we have turned on MOST BCs + * @param[in] horizontal_only don't compute vertical diffusion terms */ void DiffusionSrcForState_T (const Box& bx, const Box& domain, @@ -76,7 +77,8 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, const Array4& tm_arr, const GpuArray grav_gpu, const BCRec* bc_ptr, - const bool use_most) + const bool use_most, + const bool horizontal_only) { BL_PROFILE_VAR("DiffusionSrcForState_T()",DiffusionSrcForState_T); @@ -273,65 +275,67 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz ); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; - - Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; - rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index]) - + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) ); - - Real met_h_zeta = az(i,j,k); - - Real GradCz; - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim) ) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); - - if (ext_dir_on_zlo) { - GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + const int prim_scal_index = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? PrimScalar_comp : prim_index; + + Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_scal_index]; + rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_scal_index]) + + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_scal_index]) ); + + Real met_h_zeta = az(i,j,k); + + Real GradCz; + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == 0); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim) ) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); + + if (ext_dir_on_zlo) { + GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ); - } else if (ext_dir_on_zhi) { - GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ); - } else { - GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); - } - - if (most_on_zlo && (qty_index == RhoTheta_comp)) { - // set the exact value from MOST, don't need finite diff - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; - } + - (1./3.) * cell_prim(i, j, k+1, prim_index) ); + } else if (ext_dir_on_zhi) { + GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ); + } else { + GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); + } - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + if (most_on_zlo && (qty_index == RhoTheta_comp)) { + // set the exact value from MOST, don't need finite diff + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } - }); + }); + } // Constant rho*alpha & Turb model } else if (l_turb) { ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -394,64 +398,66 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz ); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - - Real rhoAlpha = d_alpha_eff[prim_index]; - - rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) - + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); - - Real met_h_zeta = az(i,j,k); - - Real GradCz; - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); - - if (ext_dir_on_zlo) { - GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + + Real rhoAlpha = d_alpha_eff[prim_index]; + + rhoAlpha += 0.5 * ( mu_turb(i, j, k , d_eddy_diff_idz[prim_index]) + + mu_turb(i, j, k-1, d_eddy_diff_idz[prim_index]) ); + + Real met_h_zeta = az(i,j,k); + + Real GradCz; + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == 0); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); + + if (ext_dir_on_zlo) { + GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ); - } else if (ext_dir_on_zhi) { - GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ); - } else { - GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); - } - - if (most_on_zlo && (qty_index == RhoTheta_comp)) { - // set the exact value from MOST, don't need finite diff - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; - } + - (1./3.) * cell_prim(i, j, k+1, prim_index) ); + } else if (ext_dir_on_zhi) { + GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ); + } else { + GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); + } - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + if (most_on_zlo && (qty_index == RhoTheta_comp)) { + // set the exact value from MOST, don't need finite diff + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } - }); + }); + } // Constant alpha & no LES/PBL model } else if(l_consA) { ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -512,62 +518,64 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz ); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - - Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); - Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; - - Real met_h_zeta = az(i,j,k); - - Real GradCz; - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); - - if (ext_dir_on_zlo) { - GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + + Real rhoFace = 0.5 * ( cell_data(i, j, k, Rho_comp) + cell_data(i, j, k-1, Rho_comp) ); + Real rhoAlpha = rhoFace * d_alpha_eff[prim_index]; + + Real met_h_zeta = az(i,j,k); + + Real GradCz; + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == 0); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); + + if (ext_dir_on_zlo) { + GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ); - } else if (ext_dir_on_zhi) { - GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ); - } else { - GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); - } - - if (most_on_zlo && (qty_index == RhoTheta_comp)) { - // set the exact value from MOST, don't need finite diff - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; - } + - (1./3.) * cell_prim(i, j, k+1, prim_index) ); + } else if (ext_dir_on_zhi) { + GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ); + } else { + GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); + } - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + if (most_on_zlo && (qty_index == RhoTheta_comp)) { + // set the exact value from MOST, don't need finite diff + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } - }); + }); + } // Constant rho*alpha & no LES/PBL model } else { ParallelFor(xbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept @@ -628,65 +636,69 @@ DiffusionSrcForState_T (const Box& bx, const Box& domain, yflux(i,j,k,qty_index) = -rhoAlpha * mf_v(i,j,0) * ( GradCy - (met_h_eta/met_h_zeta)*GradCz ); } }); - ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept - { - const int qty_index = start_comp + n; - const int prim_index = qty_index - 1; - - Real rhoAlpha = d_alpha_eff[prim_index]; - - Real met_h_zeta; - met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd); - - Real GradCz; - int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? - BCVars::RhoScalar_bc_comp : qty_index; - if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); - bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) - && k == 0); - bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || - (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) - && k == dom_hi.z+1); - bool most_on_zlo = ( use_most && exp_most && - (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); - - if (ext_dir_on_zlo) { - GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + if (!horizontal_only) { + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + + Real rhoAlpha = d_alpha_eff[prim_index]; + + Real met_h_zeta; + met_h_zeta = Compute_h_zeta_AtKface(i,j,k,cellSizeInv,z_nd); + + Real GradCz; + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == 0); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); + + if (ext_dir_on_zlo) { + GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + 3. * cell_prim(i, j, k , prim_index) - - (1./3.) * cell_prim(i, j, k+1, prim_index) ); - } else if (ext_dir_on_zhi) { - GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) - - 3. * cell_prim(i, j, k-1, prim_index) - + (1./3.) * cell_prim(i, j, k-2, prim_index) ); - } else { - GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); - } - - if (most_on_zlo && (qty_index == RhoTheta_comp)) { - // set the exact value from MOST, don't need finite diff - zflux(i,j,k,qty_index) = hfx_z(i,j,0); - } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { - zflux(i,j,k,qty_index) = qfx1_z(i,j,0); - } else { - zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; - } + - (1./3.) * cell_prim(i, j, k+1, prim_index) ); + } else if (ext_dir_on_zhi) { + GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ); + } else { + GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); + } - if (qty_index == RhoTheta_comp) { - if (!most_on_zlo) { - hfx_z(i,j,k) = zflux(i,j,k,qty_index); + if (most_on_zlo && (qty_index == RhoTheta_comp)) { + // set the exact value from MOST, don't need finite diff + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } else if (most_on_zlo && (qty_index == RhoQ1_comp)) { + zflux(i,j,k,qty_index) = qfx1_z(i,j,0); + } else { + zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta; } - } else if (qty_index == RhoQ1_comp) { - if (!most_on_zlo) { - qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ1_comp) { + if (!most_on_zlo) { + qfx1_z(i,j,k) = zflux(i,j,k,qty_index); + } + } else if (qty_index == RhoQ2_comp) { + qfx2_z(i,j,k) = zflux(i,j,k,qty_index); } - } else if (qty_index == RhoQ2_comp) { - qfx2_z(i,j,k) = zflux(i,j,k,qty_index); - } - }); + }); + } } // Linear combinations for z-flux with terrain + // NOTE: with horizontal_only we don't compute verticial diffusion terms, + // but do include the vertical effect from x/y diffusion with terrain //----------------------------------------------------------------------------------- // Extrapolate top and bottom cells { diff --git a/Source/ERF_IndexDefines.H b/Source/ERF_IndexDefines.H index 811946464..59d9b7a9b 100644 --- a/Source/ERF_IndexDefines.H +++ b/Source/ERF_IndexDefines.H @@ -153,7 +153,10 @@ namespace EddyDiff { Scalar_v, Q_v, PBL_lengthscale, - NumDiffs + NumDiffs, + Mom_ent_YSU = NumDiffs, + Theta_ent_YSU, + NumDiffsYSU }; } diff --git a/Source/ERF_make_new_arrays.cpp b/Source/ERF_make_new_arrays.cpp index 68630c53a..de15b5355 100644 --- a/Source/ERF_make_new_arrays.cpp +++ b/Source/ERF_make_new_arrays.cpp @@ -362,6 +362,7 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap (solverChoice.turbChoice[lev].pbl_type != PBLType::None) ); bool l_use_kturb = ( (solverChoice.turbChoice[lev].les_type != LESType::None) || (solverChoice.turbChoice[lev].pbl_type != PBLType::None) ); + bool l_use_ysu_pbl = solverChoice.turbChoice[lev].pbl_type == PBLType::YSU; bool l_use_ddorf = ( solverChoice.turbChoice[lev].les_type == LESType::Deardorff); bool l_use_moist = ( solverChoice.moisture_type != MoistureType::None ); @@ -427,7 +428,8 @@ ERF::update_diffusive_arrays (int lev, const BoxArray& ba, const DistributionMap } if (l_use_kturb) { - eddyDiffs_lev[lev] = std::make_unique(ba, dm, EddyDiff::NumDiffs, 2); + int numdiff = l_use_ysu_pbl ? EddyDiff::NumDiffsYSU : EddyDiff::NumDiffs; + eddyDiffs_lev[lev] = std::make_unique(ba, dm, numdiff, 2); eddyDiffs_lev[lev]->setVal(0.0); if(l_use_ddorf) { SmnSmn_lev[lev] = std::make_unique( ba, dm, 1, 0 ); diff --git a/Source/PBL/ERF_ComputeDiffusionSrcYSU.cpp b/Source/PBL/ERF_ComputeDiffusionSrcYSU.cpp new file mode 100644 index 000000000..4195beff7 --- /dev/null +++ b/Source/PBL/ERF_ComputeDiffusionSrcYSU.cpp @@ -0,0 +1,155 @@ +#include "ERF_ABLMost.H" +#include "ERF_DirectionSelector.H" +#include "ERF_Diffusion.H" +#include "ERF_Constants.H" +#include "ERF_TurbStruct.H" +#include "ERF_PBLModels.H" + +using namespace amrex; + +void +DiffusionSrcForStateYSU (const Box& bx, const Box& domain, + int start_comp, int num_comp, + const bool& exp_most, + const bool& rot_most, + const Array4& u, + const Array4& v, + const Array4& cell_data, + const Array4& cell_prim, + const Array4& cell_rhs, + const Array4& /*xflux*/, + const Array4& /*yflux*/, + const Array4& zflux, + const Array4& z_nd, + const Array4& /*ax*/, + const Array4& /*ay*/, + const Array4& az, + const Array4& detJ, + const GpuArray& cellSizeInv, + const Array4& SmnSmn_a, + const Array4& mf_m, + const Array4& mf_u, + const Array4& mf_v, + Array4< Real>& /*hfx_x*/, + Array4< Real>& /*hfx_y*/, + Array4< Real>& hfx_z, + Array4< Real>& /*qfx1_x*/, + Array4< Real>& /*qfx1_y*/, + Array4< Real>& /*qfx1_z*/, + Array4< Real>& /*qfx2_z*/, + Array4< Real>& diss, + const Array4& mu_turb, + const SolverChoice &solverChoice, + const int level, + const Array4& tm_arr, + const GpuArray grav_gpu, + const BCRec* bc_ptr, + const bool use_most) +{ + BL_PROFILE_VAR("DiffusionSrcForStateYSU()",DiffusionSrcForStateYSU); + + // YSU SRC: Diffusion only in the z-direction + Print() << "Computing YSU diffusion SRC for state" << std::endl; + + // For now, only theta is diffused. TODO: Moisture + int eddy_diff_idz[1] {EddyDiff::Theta_v}; + int entrainment_idz[1] {EddyDiff::Theta_ent_YSU}; + if (num_comp != 1 ) { + Abort("DiffusionSrcForStateYSU(): num_comp must be 1"); + } else if ( start_comp != RhoTheta_comp) { + Abort("DiffusionSrcForStateYSU(): start_comp must be RhoTheta_comp"); + } + + // Geometry preparations + const Box zbx = surroundingNodes(bx,2); + const auto& dom_hi = ubound(domain); + const Real dz_inv = cellSizeInv[2]; + + ParallelFor(zbx, num_comp,[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept + { + const int qty_index = start_comp + n; + const int prim_index = qty_index - 1; + const int diff_idx = eddy_diff_idz[prim_index]; + const int entrainment_idx = eddy_diff_idz[prim_index]; + + const amrex::Real rhoAlpha = 0.5 * ( mu_turb(i, j, k , diff_idx) + + mu_turb(i, j, k-1, diff_idx) ); + + const amrex::Real entrainment = 0.5 * ( mu_turb(i, j, k , entrainment_idx) + + mu_turb(i, j, k-1, entrainment_idx) ); + + Real met_h_zeta = az(i,j,k); + + Real GradCz; + int bc_comp = (qty_index >= RhoScalar_comp && qty_index < RhoScalar_comp+NSCALARS) ? + BCVars::RhoScalar_bc_comp : qty_index; + if (bc_comp > BCVars::RhoScalar_bc_comp) bc_comp -= (NSCALARS-1); + bool ext_dir_on_zlo = ( ((bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(2) == ERFBCType::ext_dir_prim)) + && k == 0); + bool ext_dir_on_zhi = ( ((bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir) || + (bc_ptr[bc_comp].lo(5) == ERFBCType::ext_dir_prim)) + && k == dom_hi.z+1); + bool most_on_zlo = ( use_most && exp_most && + (bc_ptr[bc_comp].lo(2) == ERFBCType::foextrap) && k == 0); + + if (ext_dir_on_zlo) { + GradCz = dz_inv * ( -(8./3.) * cell_prim(i, j, k-1, prim_index) + + 3. * cell_prim(i, j, k , prim_index) + - (1./3.) * cell_prim(i, j, k+1, prim_index) ); + } else if (ext_dir_on_zhi) { + GradCz = dz_inv * ( (8./3.) * cell_prim(i, j, k , prim_index) + - 3. * cell_prim(i, j, k-1, prim_index) + + (1./3.) * cell_prim(i, j, k-2, prim_index) ); + } else { + GradCz = dz_inv * ( cell_prim(i, j, k, prim_index) - cell_prim(i, j, k-1, prim_index) ); + } + + // Always True for now + if (qty_index == RhoTheta_comp) { + if (!most_on_zlo) { + zflux(i,j,k,qty_index) = -rhoAlpha * GradCz / met_h_zeta - entrainment; + hfx_z(i,j,k) = zflux(i,j,k,qty_index); + } else { + zflux(i,j,k,qty_index) = hfx_z(i,j,0); + } + } + }); + + + // Use fluxes to compute RHS + //----------------------------------------------------------------------------------- + for (int n(0); n < num_comp; ++n) + { + int qty_index = start_comp + n; + ParallelFor(bx,[=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real stateContrib = (zflux(i ,j ,k+1,qty_index) - zflux(i, j, k, qty_index)) * dz_inv; // Diffusive flux in z-dir + + stateContrib /= detJ(i,j,k); + + cell_rhs(i,j,k,qty_index) -= stateContrib; + }); + } +} + +void +DiffusionSrcForMomYSU (const Box& bxcc, const Box& bxxz , const Box& bxyz, + const Array4& uvel , + const Array4& vvel , + const Array4& wvel , + const Array4& tau33, + const Array4& tau13, const Array4& tau23, + const Array4& tau31, const Array4& tau32, + const Array4& detJ , + const GpuArray& dxInv, + const Array4& mf_m, + const Array4& /*mf_u*/, + const Array4& /*mf_v*/) +{ + BL_PROFILE_VAR("DiffusionSrcForMomYSU",DiffusionSrcForMomYSU); + Print() << "Computing YSU diffusion SRC for mom" << std::endl; + + // YSU PBL: always only vertical diffusion + // For now - only diffuse x and y velocity +} diff --git a/Source/PBL/ERF_ComputeDiffusivityYSU.cpp b/Source/PBL/ERF_ComputeDiffusivityYSU.cpp index 92503e5f7..6e74640e2 100644 --- a/Source/PBL/ERF_ComputeDiffusivityYSU.cpp +++ b/Source/PBL/ERF_ComputeDiffusivityYSU.cpp @@ -239,6 +239,11 @@ ComputeDiffusivityYSU (const MultiFab& xvel, K_turb(i,j,k,EddyDiff::Mom_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Mom_v) ,rhoKmax), rhoKmin); K_turb(i,j,k,EddyDiff::Theta_v) = std::max(std::min(K_turb(i,j,k,EddyDiff::Theta_v) ,rhoKmax), rhoKmin); K_turb(i,j,k,EddyDiff::PBL_lengthscale) = pblh_arr(i,j,0); + + // Entrainment factors: placeholder for now + K_turb(i,j,k,EddyDiff::Mom_ent_YSU) = 0.0; + K_turb(i,j,k,EddyDiff::Theta_ent_YSU) = 0.0; + }); // HACK set bottom ghost cell to 1st cell diff --git a/Source/PBL/ERF_PBLModels.H b/Source/PBL/ERF_PBLModels.H index d2fdd9d3c..3557328ab 100644 --- a/Source/PBL/ERF_PBLModels.H +++ b/Source/PBL/ERF_PBLModels.H @@ -66,6 +66,59 @@ ComputeDiffusivityYSU (const amrex::MultiFab& xvel, bool /*vert_only*/, const std::unique_ptr& z_phys_nd); +// Seperate Function to update state for YSU diffusion +void +DiffusionSrcForStateYSU (const amrex::Box& bx, const amrex::Box& domain, + int start_comp, int num_comp, + const bool& exp_most, + const bool& rot_most, + const amrex::Array4& u, + const amrex::Array4& v, + const amrex::Array4& cell_data, + const amrex::Array4& cell_prim, + const amrex::Array4& cell_rhs, + const amrex::Array4& xflux, + const amrex::Array4& yflux, + const amrex::Array4& zflux, + const amrex::Array4& z_nd, + const amrex::Array4& ax, + const amrex::Array4& ay, + const amrex::Array4& az, + const amrex::Array4& detJ, + const amrex::GpuArray& dxInv, + const amrex::Array4& SmnSmn_a, + const amrex::Array4& mf_m, + const amrex::Array4& mf_u, + const amrex::Array4& mf_v , + amrex::Array4< amrex::Real>& hfx_x, + amrex::Array4< amrex::Real>& hfx_y, + amrex::Array4< amrex::Real>& hfx_z, + amrex::Array4< amrex::Real>& qfx1_x, + amrex::Array4< amrex::Real>& qfx1_y, + amrex::Array4< amrex::Real>& qfx1_z, + amrex::Array4< amrex::Real>& qfx2_z, + amrex::Array4< amrex::Real>& diss, + const amrex::Array4& mu_turb, + const SolverChoice& solverChoice, + const int level, + const amrex::Array4& tm_arr, + const amrex::GpuArray grav_gpu, + const amrex::BCRec* bc_ptr, + const bool use_most); + +void ComputeStressesYSU (const amrex::Box& bxcc, const amrex::Box& bxxz , const amrex::Box& bxyz, + const amrex::Array4& uvel , + const amrex::Array4& vvel , + const amrex::Array4& wvel , + const amrex::Array4& tau33, + const amrex::Array4& tau13, const amrex::Array4& tau23, + const amrex::Array4& tau31, const amrex::Array4& tau32, + const amrex::Array4& detJ , + const amrex::GpuArray& dxInv, + const amrex::Array4& mf_m, + const amrex::Array4& /*mf_u*/, + const amrex::Array4& /*mf_v*/); + /** * Function for computing vertical derivatives for use in PBL model @@ -215,4 +268,5 @@ ComputeQKESourceTerms (int i, int j, int k, return source_term; } + #endif diff --git a/Source/PBL/Make.package b/Source/PBL/Make.package index 83af2fd97..92583afbc 100644 --- a/Source/PBL/Make.package +++ b/Source/PBL/Make.package @@ -1,5 +1,6 @@ CEXE_headers += ERF_PBLModels.H CEXE_sources += ERF_ComputeDiffusivityMYNN25.cpp CEXE_sources += ERF_ComputeDiffusivityYSU.cpp +CEXE_sources += ERF_ComputeDiffusionSrcYSU.cpp CEXE_headers += ERF_PBLHeight.H diff --git a/Source/TimeIntegration/ERF_TI_slow_headers.H b/Source/TimeIntegration/ERF_TI_slow_headers.H index b4959ac65..221172d73 100644 --- a/Source/TimeIntegration/ERF_TI_slow_headers.H +++ b/Source/TimeIntegration/ERF_TI_slow_headers.H @@ -10,6 +10,7 @@ #include #include +#include #include #include #include diff --git a/Source/TimeIntegration/ERF_make_tau_terms.cpp b/Source/TimeIntegration/ERF_make_tau_terms.cpp index 30c8251b0..ae02fd273 100644 --- a/Source/TimeIntegration/ERF_make_tau_terms.cpp +++ b/Source/TimeIntegration/ERF_make_tau_terms.cpp @@ -50,6 +50,7 @@ void erf_make_tau_terms (int level, int nrk, tc.les_type == LESType::Deardorff || tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); + const bool l_use_ysu_pbl = tc.pbl_type == PBLType::YSU; const bool use_most = (most != nullptr); const bool exp_most = (solverChoice.use_explicit_most); @@ -280,7 +281,7 @@ void erf_make_tau_terms (int level, int nrk, s12, s13, s21, s23, s31, s32, - er_arr, z_nd, detJ_arr, dxInv); + er_arr, z_nd, detJ_arr, dxInv, l_use_ysu_pbl); } // Remove halo cells from tau_ii but extend across valid_box bdry @@ -394,7 +395,7 @@ void erf_make_tau_terms (int level, int nrk, cell_data, s11, s22, s33, s12, s13, s23, - er_arr); + er_arr, l_use_ysu_pbl); } // Remove halo cells from tau_ii but extend across valid_box bdry diff --git a/Source/TimeIntegration/ERF_slow_rhs_post.cpp b/Source/TimeIntegration/ERF_slow_rhs_post.cpp index a91e7bdde..ebfa8ca64 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_post.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_post.cpp @@ -130,6 +130,7 @@ void erf_slow_rhs_post (int level, int finest_level, tc.les_type == LESType::Deardorff || tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); + const bool l_use_ysu_pbl = tc.pbl_type == PBLType::YSU; const bool exp_most = (solverChoice.use_explicit_most); const bool rot_most = (solverChoice.use_rotate_most); @@ -419,7 +420,7 @@ void erf_slow_rhs_post (int level, int finest_level, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z,q2fx_z, diss, mu_turb, solverChoice, level, - tm_arr, grav_gpu, bc_ptr_d, use_most); + tm_arr, grav_gpu, bc_ptr_d, use_most, l_use_ysu_pbl); } else { DiffusionSrcForState_N(tbx, domain, start_comp, num_comp, exp_most, u, v, new_cons, cur_prim, cell_rhs, @@ -427,7 +428,7 @@ void erf_slow_rhs_post (int level, int finest_level, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, - tm_arr, grav_gpu, bc_ptr_d, use_most); + tm_arr, grav_gpu, bc_ptr_d, use_most, l_use_ysu_pbl); } } // use_diff } // valid slow var diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 7ec8b392c..e9b06beb2 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -154,6 +154,7 @@ void erf_slow_rhs_pre (int level, int finest_level, tc.les_type == LESType::Deardorff || tc.pbl_type == PBLType::MYNN25 || tc.pbl_type == PBLType::YSU ); + const bool l_use_ysu_pbl = tc.pbl_type == PBLType::YSU; const bool l_use_moisture = (solverChoice.moisture_type != MoistureType::None); const bool l_use_most = (most != nullptr); @@ -507,7 +508,7 @@ void erf_slow_rhs_pre (int level, int finest_level, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, - tm_arr, grav_gpu, bc_ptr_d, l_use_most); + tm_arr, grav_gpu, bc_ptr_d, l_use_most, l_use_ysu_pbl); } else { DiffusionSrcForState_N(bx, domain, n_start, n_comp, l_exp_most, u, v, cell_data, cell_prim, cell_rhs, @@ -515,6 +516,16 @@ void erf_slow_rhs_pre (int level, int finest_level, dxInv, SmnSmn_a, mf_m, mf_u, mf_v, hfx_z, q1fx_z, q2fx_z, diss, mu_turb, solverChoice, level, + tm_arr, grav_gpu, bc_ptr_d, l_use_most, l_use_ysu_pbl); + } + if (l_use_ysu_pbl) { + DiffusionSrcForStateYSU(bx, domain, RhoTheta_comp, 1, l_exp_most, l_rot_most, u, v, + cell_data, cell_prim, cell_rhs, + diffflux_x, diffflux_y, diffflux_z, + z_nd, ax_arr, ay_arr, az_arr, detJ_arr, + dxInv, SmnSmn_a, mf_m, mf_u, mf_v, + hfx_x, hfx_y, hfx_z, q1fx_x, q1fx_y, q1fx_z, q2fx_z, diss, + mu_turb, solverChoice, level, tm_arr, grav_gpu, bc_ptr_d, l_use_most); } }