From c769d1efd69ef6f5f1a5bcee97c6d908640fd530 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Wed, 15 Nov 2023 15:29:42 -0800 Subject: [PATCH] Minor fixes to average things to their appropriate locations. (#1298) Co-authored-by: Aaron Lattanzi --- .../Diffusion/ComputeTurbulentViscosity.cpp | 4 +- Source/Diffusion/DiffusionSrcForState_N.cpp | 2 + Source/Diffusion/DiffusionSrcForState_T.cpp | 2 + Source/IO/ERF_Write1DProfiles.cpp | 9 ++-- Source/TimeIntegration/ERF_slow_rhs_pre.cpp | 42 ++++++++++--------- 5 files changed, 35 insertions(+), 24 deletions(-) diff --git a/Source/Diffusion/ComputeTurbulentViscosity.cpp b/Source/Diffusion/ComputeTurbulentViscosity.cpp index 2cc5d5946..75bc90968 100644 --- a/Source/Diffusion/ComputeTurbulentViscosity.cpp +++ b/Source/Diffusion/ComputeTurbulentViscosity.cpp @@ -126,8 +126,8 @@ void ComputeTurbulentViscosityLES (const amrex::MultiFab& Tau11, const amrex::Mu // Calculate stratification-dependent mixing length (Deardorff 1980) Real eps = std::numeric_limits::epsilon(); - Real dtheta_dz = 0.5*( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) - - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp))*dxInv[2]; + Real dtheta_dz = 0.5 * ( cell_data(i,j,k+1,RhoTheta_comp)/cell_data(i,j,k+1,Rho_comp) + - cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) )*dxInv[2]; Real E = cell_data(i,j,k,RhoKE_comp) / cell_data(i,j,k,Rho_comp); Real strat = l_abs_g * dtheta_dz * l_inv_theta0; // stratification Real length; diff --git a/Source/Diffusion/DiffusionSrcForState_N.cpp b/Source/Diffusion/DiffusionSrcForState_N.cpp index 82f884f3b..0df9470f0 100644 --- a/Source/Diffusion/DiffusionSrcForState_N.cpp +++ b/Source/Diffusion/DiffusionSrcForState_N.cpp @@ -368,6 +368,8 @@ DiffusionSrcForState_N (const amrex::Box& bx, const amrex::Box& domain, 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 + + if (qty_index==RhoTheta_comp) hfx_z(i,j,k) = -0.5 * ( zflux(i, j, k, qty_index) + zflux(i, j, k+1, qty_index) ); }); } diff --git a/Source/Diffusion/DiffusionSrcForState_T.cpp b/Source/Diffusion/DiffusionSrcForState_T.cpp index 4834d6b58..edf0a8d21 100644 --- a/Source/Diffusion/DiffusionSrcForState_T.cpp +++ b/Source/Diffusion/DiffusionSrcForState_T.cpp @@ -547,6 +547,8 @@ DiffusionSrcForState_T (const amrex::Box& bx, const amrex::Box& domain, stateContrib /= detJ(i,j,k); cell_rhs(i,j,k,qty_index) += stateContrib; + + if (qty_index==RhoTheta_comp) hfx_z(i,j,k) = -0.5 * ( zflux(i, j, k, qty_index) + zflux(i, j, k+1, qty_index) ); }); } diff --git a/Source/IO/ERF_Write1DProfiles.cpp b/Source/IO/ERF_Write1DProfiles.cpp index c34ac5937..37c5c401d 100644 --- a/Source/IO/ERF_Write1DProfiles.cpp +++ b/Source/IO/ERF_Write1DProfiles.cpp @@ -320,10 +320,13 @@ ERF::derive_stress_profiles(Gpu::HostVector& h_avg_tau11, Gpu::HostVector< ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { fab_arr(i, j, k, 0) = tau11_arr(i,j,k); - fab_arr(i, j, k, 1) = tau12_arr(i,j,k); - fab_arr(i, j, k, 2) = tau13_arr(i,j,k); + fab_arr(i, j, k, 1) = 0.25 * ( tau12_arr(i,j ,k) + tau12_arr(i+1,j ,k) + + tau12_arr(i,j+1,k) + tau12_arr(i+1,j+1,k) ); + fab_arr(i, j, k, 2) = 0.25 * ( tau13_arr(i,j,k ) + tau13_arr(i+1,j,k) + + tau13_arr(i,j,k+1) + tau13_arr(i+1,j,k+1) ); fab_arr(i, j, k, 3) = tau22_arr(i,j,k); - fab_arr(i, j, k, 4) = tau23_arr(i,j,k); + fab_arr(i, j, k, 4) = 0.25 * ( tau23_arr(i,j,k ) + tau23_arr(i,j+1,k) + + tau23_arr(i,j,k+1) + tau23_arr(i,j+1,k+1) ); fab_arr(i, j, k, 5) = tau33_arr(i,j,k); fab_arr(i, j, k, 6) = hfx3_arr(i,j,k); fab_arr(i, j, k, 7) = diss_arr(i,j,k); diff --git a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp index 08f4396d5..311532167 100644 --- a/Source/TimeIntegration/ERF_slow_rhs_pre.cpp +++ b/Source/TimeIntegration/ERF_slow_rhs_pre.cpp @@ -787,9 +787,8 @@ void erf_slow_rhs_pre (int level, int finest_level, amrex::ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation - // Add pressure gradient - amrex::Real gpx; + Real rho_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); Real met_h_xi = Compute_h_xi_AtIface (i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtIface(i, j, k, dxInv, z_nd); @@ -809,7 +808,7 @@ void erf_slow_rhs_pre (int level, int finest_level, pp_arr(i-1,j,k+1) + pp_arr(i,j,k+1) - pp_arr(i-1,j,k-1) - pp_arr(i,j,k-1) ); } - gpx = gp_xi - (met_h_xi/ met_h_zeta) * gp_zeta_on_iface; + Real gpx = gp_xi - (met_h_xi/ met_h_zeta) * gp_zeta_on_iface; gpx *= mf_u(i,j,0); Real q = 0.0; @@ -821,7 +820,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); #endif rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i-1,j,k,Rho_comp)) * abl_geo_forcing[0]; + + rho_u_face * abl_geo_forcing[0]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -834,8 +833,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_U) { - Real uu = rho_u(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_u_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (uu - dptr_rayleigh_ubar[k]) * cell_data(i,j,k,Rho_comp); + Real uu = rho_u(i,j,k) / rho_u_face; + rho_u_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (uu - dptr_rayleigh_ubar[k]) * rho_u_face; } if (l_moving_terrain) { @@ -852,6 +851,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // x-momentum equation + Real rho_u_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i-1,j,k,Rho_comp) ); Real gpx = dxInv[0] * (pp_arr(i,j,k) - pp_arr(i-1,j,k)); gpx *= mf_u(i,j,0); @@ -864,7 +864,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i-1,j,k,PrimQc_comp) ); #endif rho_u_rhs(i, j, k) += (-gpx - abl_pressure_grad[0]) / (1.0 + q) - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i-1,j,k,Rho_comp)) * abl_geo_forcing[0]; + + rho_u_face * abl_geo_forcing[0]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -895,6 +895,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation + Real rho_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); Real met_h_eta = Compute_h_eta_AtJface (i, j, k, dxInv, z_nd); Real met_h_zeta = Compute_h_zeta_AtJface(i, j, k, dxInv, z_nd); @@ -927,7 +928,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); #endif rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j-1,k,Rho_comp)) * abl_geo_forcing[1]; + + rho_v_face * abl_geo_forcing[1]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) { @@ -938,8 +939,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_V) { - Real vv = rho_v(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_v_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (vv - dptr_rayleigh_vbar[k]) * cell_data(i,j,k,Rho_comp); + Real vv = rho_v(i,j,k) / rho_v_face; + rho_v_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (vv - dptr_rayleigh_vbar[k]) * rho_v_face; } if (l_moving_terrain) { @@ -956,6 +957,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // y-momentum equation + Real rho_v_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j-1,k,Rho_comp) ); Real gpy = dxInv[1] * (pp_arr(i,j,k) - pp_arr(i,j-1,k)); gpy *= mf_v(i,j,0); @@ -968,7 +970,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j-1,k,PrimQc_comp) ); #endif rho_v_rhs(i, j, k) += (-gpy - abl_pressure_grad[1]) / (1.0_rt + q) - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j-1,k,Rho_comp)) * abl_geo_forcing[1]; + + rho_v_face * abl_geo_forcing[1]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -980,8 +982,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_V) { - Real vv = rho_v(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_v_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (vv - dptr_rayleigh_vbar[k]) * cell_data(i,j,k,Rho_comp); + Real vv = rho_v(i,j,k) / rho_v_face; + rho_v_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (vv - dptr_rayleigh_vbar[k]) * rho_v_face; } }); } // no terrain @@ -1009,6 +1011,7 @@ void erf_slow_rhs_pre (int level, int finest_level, amrex::ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation + Real rho_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); Real met_h_zeta = Compute_h_zeta_AtKface(i, j, k, dxInv, z_nd); Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ) / met_h_zeta; @@ -1021,7 +1024,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); #endif rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j,k-1,Rho_comp)) * abl_geo_forcing[2]; + + rho_w_face * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -1033,8 +1036,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_W) { - Real ww = rho_w(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_w_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (ww - dptr_rayleigh_wbar[k]) * cell_data(i,j,k,Rho_comp); + Real ww = rho_w(i,j,k) / rho_w_face; + rho_w_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (ww - dptr_rayleigh_wbar[k]) * rho_w_face; } if (l_use_terrain && l_moving_terrain) { @@ -1050,6 +1053,7 @@ void erf_slow_rhs_pre (int level, int finest_level, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // z-momentum equation + Real rho_w_face = 0.5 * ( cell_data(i,j,k,Rho_comp) + cell_data(i,j,k-1,Rho_comp) ); Real gpz = dxInv[2] * ( pp_arr(i,j,k)-pp_arr(i,j,k-1) ); Real q = 0.0; @@ -1061,7 +1065,7 @@ void erf_slow_rhs_pre (int level, int finest_level, +cell_prim(i,j,k,PrimQc_comp) + cell_prim(i,j,k-1,PrimQc_comp) ); #endif rho_w_rhs(i, j, k) += (buoyancy_fab(i,j,k) - gpz - abl_pressure_grad[2]) / (1.0_rt + q) - + 0.5*(cell_data(i,j,k,Rho_comp)+cell_data(i,j,k-1,Rho_comp)) * abl_geo_forcing[2]; + + rho_w_face * abl_geo_forcing[2]; // Add Coriolis forcing (that assumes east is +x, north is +y) if (use_coriolis) @@ -1073,8 +1077,8 @@ void erf_slow_rhs_pre (int level, int finest_level, // Add Rayleigh damping if (use_rayleigh_damping && rayleigh_damp_W) { - Real ww = rho_w(i,j,k) / cell_data(i,j,k,Rho_comp); - rho_w_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (ww - dptr_rayleigh_wbar[k]) * cell_data(i,j,k,Rho_comp); + Real ww = rho_w(i,j,k) / rho_w_face; + rho_w_rhs(i, j, k) -= dptr_rayleigh_tau[k] * (ww - dptr_rayleigh_wbar[k]) * rho_w_face; } }); } // no terrain