-
Notifications
You must be signed in to change notification settings - Fork 42
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Refactor buoyancy and add moist anelastic option.
- Loading branch information
1 parent
b2a186d
commit 6259b56
Showing
3 changed files
with
319 additions
and
160 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,229 @@ | ||
| #ifndef ERF_BUOYANCY_UTILS_H_ | ||
| #define ERF_BUOYANCY_TUILS_H_ | ||
|
|
||
| #include <ERF_EOS.H> | ||
| #include <ERF_Constants.H> | ||
|
|
||
| AMREX_GPU_HOST | ||
| AMREX_FORCE_INLINE | ||
| amrex::Real | ||
| buoyancy_dry_default (int& i, | ||
| int& j, | ||
| int& k, | ||
| amrex::Real const& grav_gpu, | ||
| amrex::Real const& rd_over_cp, | ||
| const amrex::Array4<const amrex::Real>& p0_arr, | ||
| const amrex::Array4<const amrex::Real>& r0_arr, | ||
| const amrex::Array4<const amrex::Real>& cell_data) | ||
| { | ||
| amrex::Real rt0_hi = getRhoThetagivenP(p0_arr(i,j,k)); | ||
| amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), rt0_hi/r0_arr(i,j,k), rd_over_cp); | ||
| amrex::Real t_hi = getTgivenPandTh(p0_arr(i,j,k), cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k), rd_over_cp); | ||
| amrex::Real qplus = (t_hi-t0_hi)/t0_hi; | ||
|
|
||
| amrex::Real rt0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)); | ||
| amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), rt0_lo/r0_arr(i,j,k-1), rd_over_cp); | ||
| amrex::Real t_lo = getTgivenPandTh(p0_arr(i,j,k-1), cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1), rd_over_cp); | ||
| amrex::Real qminus = (t_lo-t0_lo)/t0_lo; | ||
|
|
||
| amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * qplus + r0_arr(i,j,k-1) * qminus); | ||
| return (-r0_q_avg * grav_gpu); | ||
| } | ||
|
|
||
|
|
||
| AMREX_GPU_HOST | ||
| AMREX_FORCE_INLINE | ||
| amrex::Real | ||
| buoyancy_moist_anelastic (int& i, | ||
| int& j, | ||
| int& k, | ||
| const int& klo, | ||
| amrex::Real const& grav_gpu, | ||
| amrex::Real const& rv_over_rd, | ||
| const amrex::Array4<const amrex::Real>& p0_arr, | ||
| const amrex::Array4<const amrex::Real>& r0_arr, | ||
| const amrex::Array4<const amrex::Real>& cell_data) | ||
| { | ||
| amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1); | ||
| amrex::Real qv_lo = cell_data(i,j,k-1,RhoQ1_comp)/r0_arr(i,j,k-1); | ||
| amrex::Real qc_lo = cell_data(i,j,k-1,RhoQ2_comp)/r0_arr(i,j,k-1); | ||
| amrex::Real qt_lo = qv_lo + qc_lo; | ||
| amrex::Real theta_v_lo = theta_d_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo); | ||
|
|
||
| amrex::Real theta_d_hi = cell_data(i,j,k,RhoTheta_comp)/r0_arr(i,j,k); | ||
| amrex::Real qv_hi = cell_data(i,j,k,RhoQ1_comp)/r0_arr(i,j,k); | ||
| amrex::Real qc_hi = cell_data(i,j,k,RhoQ2_comp)/r0_arr(i,j,k); | ||
| amrex::Real qt_hi = qv_hi + qc_hi; | ||
| amrex::Real theta_v_hi = theta_d_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi); | ||
|
|
||
| amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi); | ||
|
|
||
| amrex::Real theta_d_0_lo = getRhoThetagivenP(p0_arr(i,j,k-1)) / r0_arr(i,j,k-1); | ||
| amrex::Real theta_d_0_hi = getRhoThetagivenP(p0_arr(i,j,k )) / r0_arr(i,j,k ); | ||
| amrex::Real theta_v_0_lo = theta_d_0_lo * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo); | ||
| amrex::Real theta_v_0_hi = theta_d_0_hi * (1.0 - (1.0 - rv_over_rd)*qt_hi - rv_over_rd*qc_hi); | ||
|
|
||
| amrex::Real theta_v_0_wface = amrex::Real(0.5) * (theta_v_0_lo + theta_v_0_hi); | ||
|
|
||
| return (-grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface); | ||
| } | ||
|
|
||
|
|
||
| AMREX_GPU_HOST | ||
| AMREX_FORCE_INLINE | ||
| amrex::Real | ||
| buoyancy_type1 (int& i, | ||
| int& j, | ||
| int& k, | ||
| const int& n_qstate, | ||
| amrex::Real const& grav_gpu, | ||
| const amrex::Array4<const amrex::Real>& r0_arr, | ||
| const amrex::Array4<const amrex::Real>& cell_data) | ||
| { | ||
| amrex::Real rhop_hi = cell_data(i,j,k ,Rho_comp) - r0_arr(i,j,k ); | ||
| amrex::Real rhop_lo = cell_data(i,j,k-1,Rho_comp) - r0_arr(i,j,k-1); | ||
| for (int q_offset(0); q_offset<n_qstate; ++q_offset) { | ||
| rhop_hi += cell_data(i,j,k ,RhoQ1_comp+q_offset); | ||
| rhop_lo += cell_data(i,j,k-1,RhoQ1_comp+q_offset); | ||
| } | ||
| return( grav_gpu * amrex::Real(0.5) * ( rhop_hi + rhop_lo ) ); | ||
| } | ||
|
|
||
|
|
||
| AMREX_GPU_HOST | ||
| AMREX_FORCE_INLINE | ||
| amrex::Real | ||
| buoyancy_type2 (int& i, | ||
| int& j, | ||
| int& k, | ||
| const int& n_qstate, | ||
| amrex::Real const& grav_gpu, | ||
| amrex::Real* rho_d_ptr, | ||
| amrex::Real* theta_d_ptr, | ||
| amrex::Real* qv_d_ptr, | ||
| amrex::Real* qc_d_ptr, | ||
| amrex::Real* qp_d_ptr, | ||
| const amrex::Array4<const amrex::Real>& cell_prim, | ||
| const amrex::Array4<const amrex::Real>& cell_data) | ||
| { | ||
| amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]); | ||
| amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]); | ||
|
|
||
| amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), | ||
| cell_data(i,j,k ,RhoTheta_comp), | ||
| cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp)); | ||
| amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), | ||
| cell_data(i,j,k-1,RhoTheta_comp), | ||
| cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp)); | ||
|
|
||
| amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; | ||
| amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; | ||
|
|
||
| amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; | ||
| amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; | ||
|
|
||
| amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; | ||
| amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; | ||
|
|
||
| amrex::Real qplus = 0.61 * ( qv_plus - qv_d_ptr[k] ) - | ||
| ( qc_plus - qc_d_ptr[k] + | ||
| qp_plus - qp_d_ptr[k] ) | ||
| + (tempp3d-tempp1d)/tempp1d*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k]-qc_d_ptr[k]-qp_d_ptr[k]); | ||
|
|
||
| amrex::Real qminus = 0.61 * ( qv_minus - qv_d_ptr[k-1] ) - | ||
| ( qc_minus - qc_d_ptr[k-1] + | ||
| qp_minus - qp_d_ptr[k-1] ) | ||
| + (tempm3d-tempm1d)/tempm1d*(amrex::Real(1.0) + amrex::Real(0.61)*qv_d_ptr[k-1]-qc_d_ptr[k-1]-qp_d_ptr[k-1]); | ||
|
|
||
| amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus); | ||
| amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); | ||
|
|
||
| return (-qavg * r0avg * grav_gpu); | ||
| } | ||
|
|
||
|
|
||
| AMREX_GPU_HOST | ||
| AMREX_FORCE_INLINE | ||
| amrex::Real | ||
| buoyancy_type3 (int& i, | ||
| int& j, | ||
| int& k, | ||
| const int& n_qstate, | ||
| amrex::Real const& grav_gpu, | ||
| amrex::Real* rho_d_ptr, | ||
| amrex::Real* theta_d_ptr, | ||
| amrex::Real* qv_d_ptr, | ||
| const amrex::Array4<const amrex::Real>& cell_prim, | ||
| const amrex::Array4<const amrex::Real>& cell_data) | ||
| { | ||
| amrex::Real tempp1d = getTgivenRandRTh(rho_d_ptr[k ], rho_d_ptr[k ]*theta_d_ptr[k ], qv_d_ptr[k ]); | ||
| amrex::Real tempm1d = getTgivenRandRTh(rho_d_ptr[k-1], rho_d_ptr[k-1]*theta_d_ptr[k-1], qv_d_ptr[k-1]); | ||
|
|
||
| amrex::Real tempp3d = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), | ||
| cell_data(i,j,k ,RhoTheta_comp), | ||
| cell_data(i,j,k ,RhoQ1_comp)/cell_data(i,j,k ,Rho_comp)); | ||
| amrex::Real tempm3d = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), | ||
| cell_data(i,j,k-1,RhoTheta_comp), | ||
| cell_data(i,j,k-1,RhoQ1_comp)/cell_data(i,j,k-1,Rho_comp)); | ||
|
|
||
| amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; | ||
| amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; | ||
|
|
||
| amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; | ||
| amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; | ||
|
|
||
| amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; | ||
| amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; | ||
|
|
||
| amrex::Real qplus = 0.61 * qv_plus - (qc_plus + qp_plus) + (tempp3d-tempp1d)/tempp1d; | ||
|
|
||
| amrex::Real qminus = 0.61 * qv_minus - (qc_minus + qp_minus) + (tempm3d-tempm1d)/tempm1d; | ||
|
|
||
| amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus); | ||
| amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); | ||
|
|
||
| return ( -qavg * r0avg * grav_gpu ); | ||
| } | ||
|
|
||
|
|
||
| AMREX_GPU_HOST | ||
| AMREX_FORCE_INLINE | ||
| amrex::Real | ||
| buoyancy_type4 (int& i, | ||
| int& j, | ||
| int& k, | ||
| const int& n_qstate, | ||
| amrex::Real const& grav_gpu, | ||
| amrex::Real* rho_d_ptr, | ||
| amrex::Real* theta_d_ptr, | ||
| amrex::Real* qv_d_ptr, | ||
| amrex::Real* qc_d_ptr, | ||
| amrex::Real* qp_d_ptr, | ||
| const amrex::Array4<const amrex::Real>& cell_prim, | ||
| const amrex::Array4<const amrex::Real>& cell_data) | ||
| { | ||
| amrex::Real qv_plus = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; | ||
| amrex::Real qv_minus = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; | ||
|
|
||
| amrex::Real qc_plus = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; | ||
| amrex::Real qc_minus = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; | ||
|
|
||
| amrex::Real qp_plus = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; | ||
| amrex::Real qp_minus = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_comp) : 0.0; | ||
|
|
||
| amrex::Real qplus = amrex::Real(0.61) * ( qv_plus - qv_d_ptr[k] ) - | ||
| ( qc_plus - qc_d_ptr[k] + | ||
| qp_plus - qp_d_ptr[k] ) | ||
| + (cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp) - theta_d_ptr[k ])/theta_d_ptr[k ]; | ||
|
|
||
| amrex::Real qminus = amrex::Real(0.61) * ( qv_minus - qv_d_ptr[k-1] ) - | ||
| ( qc_minus - qc_d_ptr[k-1] + | ||
| qp_minus - qp_d_ptr[k-1] ) | ||
| + (cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp) - theta_d_ptr[k-1])/theta_d_ptr[k-1]; | ||
|
|
||
| amrex::Real qavg = amrex::Real(0.5) * (qplus + qminus); | ||
| amrex::Real r0avg = amrex::Real(0.5) * (rho_d_ptr[k] + rho_d_ptr[k-1]); | ||
|
|
||
| return (-qavg * r0avg * grav_gpu); | ||
| } | ||
| #endif |
Oops, something went wrong.