Skip to content

Commit

Permalink
Subsidence w/ terrain (erf-model#1944)
Browse files Browse the repository at this point in the history 
* Account for grid stretching with subsidence term.

* Handle dzInv more rigorously with terrain.

---------

Co-authored-by: Aaron Lattanzi <[email protected]>
  • Loading branch information
@AMLattanzi
AMLattanzi and Aaron Lattanzi authored Nov 13, 2024
1 parent 50fc1a5 commit 75c4749
Show file tree
Hide file tree
Showing 2 changed files with 46 additions and 20 deletions.
44 changes: 36 additions & 8 deletions Source/SourceTerms/ERF_make_mom_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -342,38 +342,66 @@ void make_mom_sources (int level,
const int nr = Rho_comp;
ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = 0.5*dxInv[2];
if (z_nd_arr) {
Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
+ z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
+ z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
dzInv = 1.0 / (z_xf_hi - z_xf_lo);
}
Real rho_on_u_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i-1,j,k,nr) );
Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf *
0.5 * (U_hi - U_lo) * dxInv[2];
xmom_src_arr(i, j, k) -= rho_on_u_face * wbar_xf * (U_hi - U_lo) * dzInv;
});
ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = 0.5*dxInv[2];
if (z_nd_arr) {
Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
+ z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
+ z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
dzInv = 1.0 / (z_yf_hi - z_yf_lo);
}
Real rho_on_v_face = 0.5 * ( cell_data(i,j,k,nr) + cell_data(i,j-1,k,nr) );
Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf *
0.5 * (V_hi - V_lo) * dxInv[2];
ymom_src_arr(i, j, k) -= rho_on_v_face * wbar_yf * (V_hi - V_lo) * dzInv;
});
} else {
ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = 0.5*dxInv[2];
if (z_nd_arr) {
Real z_xf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i,j+1,k )
+ z_nd_arr(i,j,k-1) + z_nd_arr(i,j+1,k-1) );
Real z_xf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i,j+1,k+1)
+ z_nd_arr(i,j,k+2) + z_nd_arr(i,j+1,k+2) );
dzInv = 1.0 / (z_xf_hi - z_xf_lo);
}
Real U_hi = dptr_u_plane(k+1) / dptr_r_plane(k+1);
Real U_lo = dptr_u_plane(k-1) / dptr_r_plane(k-1);
Real wbar_xf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
xmom_src_arr(i, j, k) -= wbar_xf *
0.5 * (U_hi - U_lo) * dxInv[2];
xmom_src_arr(i, j, k) -= wbar_xf * (U_hi - U_lo) * dzInv;
});
ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = 0.5*dxInv[2];
if (z_nd_arr) {
Real z_yf_lo = 0.25 * ( z_nd_arr(i,j,k ) + z_nd_arr(i+1,j,k )
+ z_nd_arr(i,j,k-1) + z_nd_arr(i+1,j,k-1) );
Real z_yf_hi = 0.25 * ( z_nd_arr(i,j,k+1) + z_nd_arr(i+1,j,k+1)
+ z_nd_arr(i,j,k+2) + z_nd_arr(i+1,j,k+2) );
dzInv = 1.0 / (z_yf_hi - z_yf_lo);
}
Real V_hi = dptr_v_plane(k+1) / dptr_r_plane(k+1);
Real V_lo = dptr_v_plane(k-1) / dptr_r_plane(k-1);
Real wbar_yf = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
ymom_src_arr(i, j, k) -= wbar_yf *
0.5 * (V_hi - V_lo) * dxInv[2];
ymom_src_arr(i, j, k) -= wbar_yf * (V_hi - V_lo) * dzInv;
});
}
}
Expand Down
22 changes: 10 additions & 12 deletions Source/SourceTerms/ERF_make_sources.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -274,20 +274,20 @@ void make_sources (int level,
const int nr = Rho_comp;
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = (z_cc_arr) ? 1.0/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : 0.5*dxInv[2];
Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, n) -= cell_data(i,j,k,nr) * wbar_cc *
0.5 * (T_hi - T_lo) * dxInv[2];
cell_src(i, j, k, n) -= cell_data(i,j,k,nr) * wbar_cc * (T_hi - T_lo) * dzInv;
});
} else {
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = (z_cc_arr) ? 1.0/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : 0.5*dxInv[2];
Real T_hi = dptr_t_plane(k+1) / dptr_r_plane(k+1);
Real T_lo = dptr_t_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, n) -= wbar_cc *
0.5 * (T_hi - T_lo) * dxInv[2];
cell_src(i, j, k, n) -= wbar_cc * (T_hi - T_lo) * dzInv;
});
}
}
Expand All @@ -301,28 +301,26 @@ void make_sources (int level,
const int nr = Rho_comp;
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = (z_cc_arr) ? 1.0/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : 0.5*dxInv[2];
Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
Real Qc_hi = dptr_qc_plane(k+1) / dptr_r_plane(k+1);
Real Qc_lo = dptr_qc_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, nv ) -= cell_data(i,j,k,nr) * wbar_cc *
0.5 * (Qv_hi - Qv_lo) * dxInv[2];
cell_src(i, j, k, nv+1) -= cell_data(i,j,k,nr) * wbar_cc *
0.5 * (Qc_hi - Qc_lo) * dxInv[2];
cell_src(i, j, k, nv ) -= cell_data(i,j,k,nr) * wbar_cc * (Qv_hi - Qv_lo) * dzInv;
cell_src(i, j, k, nv+1) -= cell_data(i,j,k,nr) * wbar_cc * (Qc_hi - Qc_lo) * dzInv;
});
} else {
ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept
{
Real dzInv = (z_cc_arr) ? 1.0/ (z_cc_arr(i,j,k+1) - z_cc_arr(i,j,k-1)) : 0.5*dxInv[2];
Real Qv_hi = dptr_qv_plane(k+1) / dptr_r_plane(k+1);
Real Qv_lo = dptr_qv_plane(k-1) / dptr_r_plane(k-1);
Real Qc_hi = dptr_qc_plane(k+1) / dptr_r_plane(k+1);
Real Qc_lo = dptr_qc_plane(k-1) / dptr_r_plane(k-1);
Real wbar_cc = 0.5 * (dptr_wbar_sub[k] + dptr_wbar_sub[k+1]);
cell_src(i, j, k, nv ) -= wbar_cc *
0.5 * (Qv_hi - Qv_lo) * dxInv[2];
cell_src(i, j, k, nv+1) -= wbar_cc *
0.5 * (Qc_hi - Qc_lo) * dxInv[2];
cell_src(i, j, k, nv ) -= wbar_cc * (Qv_hi - Qv_lo) * dzInv;
cell_src(i, j, k, nv+1) -= wbar_cc * (Qc_hi - Qc_lo) * dzInv;
});
}
}
Expand Down

0 comments on commit 75c4749

Please sign in to comment.