Skip to content

Commit

Permalink
Minor fixes to average things to their appropriate locations. (#1298)
Browse files Browse the repository at this point in the history 
Co-authored-by: Aaron Lattanzi <[email protected]>
  • Loading branch information
@AMLattanzi
AMLattanzi and Aaron Lattanzi authored Nov 15, 2023
1 parent 89f00bb commit c769d1e
Show file tree
Hide file tree
Showing 5 changed files with 35 additions and 24 deletions.
4 changes: 2 additions & 2 deletions Source/Diffusion/ComputeTurbulentViscosity.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<Real>::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;
Expand Down
2 changes: 2 additions & 0 deletions Source/Diffusion/DiffusionSrcForState_N.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) );
});
}

Expand Down
2 changes: 2 additions & 0 deletions Source/Diffusion/DiffusionSrcForState_T.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) );
});
}

Expand Down
9 changes: 6 additions & 3 deletions Source/IO/ERF_Write1DProfiles.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,10 +320,13 @@ ERF::derive_stress_profiles(Gpu::HostVector<Real>& 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);
Expand Down
42 changes: 23 additions & 19 deletions Source/TimeIntegration/ERF_slow_rhs_pre.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand All @@ -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;
Expand All @@ -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)
Expand All @@ -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) {
Expand All @@ -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);

Expand All @@ -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)
Expand Down Expand Up @@ -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);

Expand Down Expand Up @@ -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)
{
Expand All @@ -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) {
Expand All @@ -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);

Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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;

Expand All @@ -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)
Expand All @@ -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) {
Expand All @@ -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;
Expand All @@ -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)
Expand All @@ -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
Expand Down

0 comments on commit c769d1e

Please sign in to comment.