diff --git a/Docs/sphinx_doc/LinearSolvers.rst b/Docs/sphinx_doc/LinearSolvers.rst index 6b6a01f2a..1d34c2bb8 100644 --- a/Docs/sphinx_doc/LinearSolvers.rst +++ b/Docs/sphinx_doc/LinearSolvers.rst @@ -17,8 +17,9 @@ Multigrid can also be used when the union of grids at a level is not in itself r For simulations using grid stretching in the vertical but flat terrain, we must use the hybrid FFT solver in which we perform 2D transforms only in the lateral directions and couple the solution in the vertical direction with a tridiagonal solve. In both these cases we use a 7-point stencil. + To solve the Poisson equation on terrain-fitted coordinates with general terrain, we rely on the FFT-preconditioned GMRES solver since the stencil effectively has variable coefficients and requires 19 points. - .. note:: - **Currently only doubly periodic lateral boundary conditions are supported by the hybrid FFT, and therefore by the GMRES solver. More general boundary conditions are a work in progress.** +- .. note:: +- **The FFT solver / preconditioner can only be used when the union of grids at a level is itself rectangular. diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index b353fa291..4d4b28ea7 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -7,11 +7,24 @@ .. _sec:Buoyancy: +ERF has several options for how to define the buoyancy force. + Buoyancy ========= -ERF has three options for how to define the buoyancy force. Even in the absence of moisture these -expressions are not equivalent. +When using the anelastic formulation, the expression for buoyancy used by ERF is + +.. math:: + \mathbf{B} = -\rho^\prime \mathbf{g} \approx -\rho_0 \mathbf{g} ( \frac{\theta_v^\prime}{\overline{\theta_0}} + +where we define + +.. math:: + \theta_v = \theta_d (1 - (1 - \frac{R_v}{R_d}) (q_v + q_c) - \frac{R_v}{R_d} q_c) + + +When using the compressible formulation, the following choices are available. + Density of the mixture ----------------------- @@ -52,56 +65,10 @@ for the background state. For eg., a usual scenario is that of a background stat the total background density :math:`\rho_0 = \rho_{d_0}(1 + q_{v_0})`, where :math:`\rho_{d_0}` and :math:`q_{v_0}` are the background dry density and vapor mixing ratio respectively. As a check, we observe that :math:`\rho^\prime_0 = 0`, which means that the background state is not buoyant. -Type 2 ------- - -The second option for the buoyancy force is - -.. math:: - \mathbf{B} = -\rho_0 \mathbf{g} ( 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime - + \frac{T^\prime}{\bar{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ) - -To derive this expression, we define :math:`T_v = T (1 + 0.61 q_v − q_c − q_i - q_p)`, then we can write - -.. math:: - p = \rho (R_d q_d + R_v q_v) T = \rho R_d T (1 + 0.61 q_v − q_c − q_i - q_p ) = \rho R_d T_v - - -Starting from :math:`p = \rho R_d T_v` and neglecting :math:`\frac{p^\prime}{\bar{p}}`, we now write - -.. math:: - \frac{\rho^\prime}{\overline{\rho}} = -\frac{T_v^\prime}{\overline{T_v}} - -and define - -.. math:: - - T_v^\prime = T_v - \overline{T_v} \approx \overline{T} [ 0.61 q_v^\prime - (q_c^\prime + q_i^\prime + q_p^\prime)] + - (T - \overline{T}) [1+ 0.61 \bar{q_v} - \bar{q_c} - \bar{q_i} - \bar{q_p} ] . - -where we have retained only first order terms in perturbational quantities. - -Then - -.. math:: - - \mathbf{B} = \rho^\prime \mathbf{g} = -\overline{\rho} \frac{\overline{T}}{\overline{T_v}} \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime - q_i^\prime - q_p^\prime ) + \frac{T^\prime}{\overline{T_v}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ] - -where the overbar represents a horizontal average of the current state and the perturbation is defined relative to that average. - -Again keeping only the first order terms in the mass mixing ratios, we can simplify this to - -.. math:: - \mathbf{B} = \rho^\prime \mathbf{g} = -\rho_0 \mathbf{g} [ 0.61 q_v^\prime - q_c^\prime + q_i^\prime + q_p^\prime - + \frac{T^\prime}{\overline{T}} (1.0 + 0.61 \bar{q_v} - \bar{q_i} - \bar{q_c} - \bar{q_p}) ] - -We note that this reduces to Type 3 if the horizontal averages of the moisture terms are all zero. - Type 3 ------ -The third formulation of the buoyancy term assumes that the horizontal averages of the moisture quantities are negligible, -which removes the need to compute horizontal averages of these quantities. This reduces the Type 2 expression to the following: +This formulation of the buoyancy term assumes that the horizontal averages of the moisture quantities are negligible. .. math:: \mathbf{B} = \rho^\prime \mathbf{g} \approx -\rho_0 \mathbf{g} ( \frac{T^\prime}{\overline{T}} @@ -120,7 +87,7 @@ This expression for buoyancy is from `khairoutdinov2003cloud`_ and `bryan2002ben .. math:: \begin{equation} - \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg)\mathbf{g}, + \mathbf{B} = \rho'\mathbf{g} \approx -\rho \mathbf{g} \Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg) \end{equation} The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as diff --git a/Source/DataStructs/ERF_DataStruct.H b/Source/DataStructs/ERF_DataStruct.H index d8f05797b..41707c247 100644 --- a/Source/DataStructs/ERF_DataStruct.H +++ b/Source/DataStructs/ERF_DataStruct.H @@ -126,14 +126,15 @@ struct SolverChoice { if (moisture_type == MoistureType::Kessler_NoRain || moisture_type == MoistureType::SAM || moisture_type == MoistureType::SAM_NoIce || - moisture_type == MoistureType::SAM_NoPrecip_NoIce) { - buoyancy_type = 1; // asserted in make buoyancy + moisture_type == MoistureType::SAM_NoPrecip_NoIce) + { + buoyancy_type = 1; // uses Rhoprime } else { buoyancy_type = 2; // uses Tprime } } - // Which expression (1,2 or 3) to use for buoyancy + // Which expression (1,2/3 or 4) to use for buoyancy pp.query("buoyancy_type", buoyancy_type); // What type of land surface model to use @@ -196,7 +197,7 @@ struct SolverChoice { if (any_anelastic == 1) { constant_density = true; project_initial_velocity = true; - buoyancy_type = 2; // uses Tprime + buoyancy_type = 3; // (This isn't actually used when anelastic is set) } else { pp.query("project_initial_velocity", project_initial_velocity); diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H index fd18827f3..967a954fc 100644 --- a/Source/SourceTerms/ERF_buoyancy_utils.H +++ b/Source/SourceTerms/ERF_buoyancy_utils.H @@ -17,16 +17,11 @@ buoyancy_dry_anelastic (int& i, int& j, int& k, amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/cell_data(i,j,k-1,Rho_comp); amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/cell_data(i,j,k ,Rho_comp); - amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi); + amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_hi); + amrex::Real theta_d0_wface = amrex::Real(0.5) * (th0_arr(i,j,k) + th0_arr(i,j,k-1)); + amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - amrex::Real theta_d_0_lo = th0_arr(i,j,k-1); - amrex::Real theta_d_0_hi = th0_arr(i,j,k ); - - amrex::Real theta_d_0_wface = amrex::Real(0.5) * (theta_d_0_lo + theta_d_0_hi); - - amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - - return (-rho0_wface * grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface); + return (-rho0_wface * grav_gpu * (theta_d_wface - theta_d0_wface) / theta_d0_wface); } AMREX_GPU_DEVICE @@ -51,39 +46,11 @@ buoyancy_moist_anelastic (int& i, int& j, int& 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_v_0_lo = th0_arr(i,j,k-1) * (1.0 - (1.0 - rv_over_rd)*qt_lo - rv_over_rd*qc_lo); - amrex::Real theta_v_0_hi = th0_arr(i,j,k ) * (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); - - amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - - return (-rho0_wface * grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface); -} - -AMREX_GPU_DEVICE -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& r0_arr, - const amrex::Array4& p0_arr, - const amrex::Array4& th0_arr, - const amrex::Array4& cell_data) -{ - amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); - amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); - amrex::Real qplus = (t_hi-t0_hi)/t0_hi; - - amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp); - amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - amrex::Real qminus = (t_lo-t0_lo)/t0_lo; + amrex::Real theta_v_wface = amrex::Real(0.5) * (theta_v_lo + theta_v_hi); + amrex::Real theta_v0_wface = amrex::Real(0.5) * (th0_arr(i,j,k) + th0_arr(i,j,k-1)); + amrex::Real rho0_wface = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - 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); + return (-rho0_wface * grav_gpu * (theta_v_wface - theta_v0_wface) / theta_v0_wface); } AMREX_GPU_DEVICE @@ -104,165 +71,138 @@ buoyancy_rhopert (int& i, int& j, int& k, return( grav_gpu * amrex::Real(0.5) * ( rhop_hi + rhop_lo ) ); } - AMREX_GPU_DEVICE 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& cell_prim, - const amrex::Array4& cell_data) +buoyancy_dry_Tpert (int& i, int& j, int& k, + amrex::Real const& grav_gpu, + amrex::Real const& rd_over_cp, + const amrex::Array4& r0_arr, + const amrex::Array4& p0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& 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 t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); + amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_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 t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp); + amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - 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 tprime_hi = (t_hi-t0_hi)/t0_hi; + amrex::Real tprime_lo = (t_lo-t0_lo)/t0_lo; - 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 tp_avg = amrex::Real(0.5) * (tprime_hi + tprime_lo); + amrex::Real r0_avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - 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]); + return ( -r0_avg * grav_gpu * tp_avg); +} - 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_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +buoyancy_dry_Thpert (int& i, int& j, int& k, + amrex::Real const& grav_gpu, + const amrex::Array4& r0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& cell_prim) +{ + // + // TODO: we are currently using Theta_prime/Theta_0 to replace T_prime/T_0 - P_prime/P_0 + // but I don't think that is quite right. + // + amrex::Real thetaprime_hi = (cell_prim(i,j,k ,PrimTheta_comp) - th0_arr(i,j,k )) / th0_arr(i,j,k ); + amrex::Real thetaprime_lo = (cell_prim(i,j,k-1,PrimTheta_comp) - th0_arr(i,j,k-1)) / th0_arr(i,j,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]); + amrex::Real thp_avg = amrex::Real(0.5) * (thetaprime_hi + thetaprime_lo); + amrex::Real r0_avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - return (-qavg * r0avg * grav_gpu); + return ( -r0_avg * grav_gpu * thp_avg); } - AMREX_GPU_DEVICE 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& cell_prim, - const amrex::Array4& cell_data) +buoyancy_moist_Tpert (int& i, int& j, int& k, + const int& n_qstate, + amrex::Real const& grav_gpu, + amrex::Real const& rd_over_cp, + const amrex::Array4& r0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& p0_arr, + const amrex::Array4& cell_prim, + const amrex::Array4& 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]); + // + // Note: this currently assumes the base state qv0 is identically zero + // TODO: generalize this to allow for moist base state + // + amrex::Real qv_hi = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; + amrex::Real qv_lo = (n_qstate >= 1) ? cell_prim(i,j,k-1,PrimQ1_comp) : 0.0; - 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 qc_hi = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; + amrex::Real qc_lo = (n_qstate >= 2) ? cell_prim(i,j,k-1,PrimQ2_comp) : 0.0; - 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 qp_hi = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; + amrex::Real qp_lo = (n_qstate >= 3) ? cell_prim(i,j,k-1,PrimQ3_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 t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); + amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); - 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 t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp), qv_hi); + amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp), qv_lo); - amrex::Real qplus = 0.61 * qv_plus - (qc_plus + qp_plus) + (tempp3d-tempp1d)/tempp1d; + amrex::Real q_hi = 0.61 * qv_hi - (qc_hi + qp_hi) + (t_hi-t0_hi)/t0_hi; + amrex::Real q_lo = 0.61 * qv_lo - (qc_lo + qp_lo) + (t_lo-t0_lo)/t0_lo; - amrex::Real qminus = 0.61 * qv_minus - (qc_minus + qp_minus) + (tempm3d-tempm1d)/tempm1d; + amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo); + amrex::Real r0avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,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 ); + return ( -r0avg * grav_gpu * qavg); } - AMREX_GPU_DEVICE 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& cell_prim, - const amrex::Array4& cell_data) +buoyancy_moist_Thpert (int& i, int& j, int& k, + const int& n_qstate, + amrex::Real const& grav_gpu, + const amrex::Array4& r0_arr, + const amrex::Array4& th0_arr, + const amrex::Array4& cell_prim) { - 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; + // + // Note: this currently assumes the base state qv0 is identically zero + // TODO: generalize this to allow for moist base state + // + amrex::Real qv_hi = (n_qstate >= 1) ? cell_prim(i,j,k ,PrimQ1_comp) : 0.0; + amrex::Real qv_lo = (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 qc_hi = (n_qstate >= 2) ? cell_prim(i,j,k ,PrimQ2_comp) : 0.0; + amrex::Real qc_lo = (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 qp_hi = (n_qstate >= 3) ? cell_prim(i,j,k ,PrimQ3_comp) : 0.0; + amrex::Real qp_lo = (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 ]; + // + // TODO: we are currently using Theta_prime/Theta_0 to replace T_prime/T_0 - P_prime/P_0 + // but I don't think that is quite right. + // + amrex::Real q_hi = amrex::Real(0.61) * qv_hi - (qc_hi + qp_hi) + + (cell_prim(i,j,k ,PrimTheta_comp) - th0_arr(i,j,k )) / th0_arr(i,j,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 q_lo = amrex::Real(0.61) * qv_lo - (qc_lo + qp_lo) + + (cell_prim(i,j,k-1,PrimTheta_comp) - th0_arr(i,j,k-1)) / th0_arr(i,j,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]); + amrex::Real qavg = amrex::Real(0.5) * (q_hi + q_lo); + amrex::Real r0avg = amrex::Real(0.5) * (r0_arr(i,j,k) + r0_arr(i,j,k-1)); - return (-qavg * r0avg * grav_gpu); + return ( -r0avg * grav_gpu * qavg); } // ************************************************************************************** // Routines below this line are not currently used // ************************************************************************************** -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -amrex::Real -buoyancy_dry_Tpert (int& i, int& j, int& k, - amrex::Real const& grav_gpu, - amrex::Real const& rd_over_cp, - const amrex::Array4& r0_arr, - const amrex::Array4& p0_arr, - const amrex::Array4& th0_arr, - const amrex::Array4& cell_data) -{ - amrex::Real t0_hi = getTgivenPandTh(p0_arr(i,j,k), th0_arr(i,j,k), rd_over_cp); - amrex::Real t_hi = getTgivenRandRTh(cell_data(i,j,k ,Rho_comp), cell_data(i,j,k ,RhoTheta_comp)); - amrex::Real qplus = (t_hi-t0_hi)/t0_hi; - - amrex::Real t0_lo = getTgivenPandTh(p0_arr(i,j,k-1), th0_arr(i,j,k-1), rd_over_cp); - amrex::Real t_lo = getTgivenRandRTh(cell_data(i,j,k-1,Rho_comp), cell_data(i,j,k-1,RhoTheta_comp)); - 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_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -276,14 +216,14 @@ buoyancy_dry_anelastic_T (int& i, int& j, int& k, 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 q_hi = (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 q_lo = (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); + amrex::Real r0_q_avg = amrex::Real(0.5) * (r0_arr(i,j,k) * q_hi + r0_arr(i,j,k-1) * q_lo); return (-r0_q_avg * grav_gpu); } diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index 416c79a71..c42882319 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -53,267 +53,133 @@ void make_buoyancy (Vector& S_data, MultiFab p0 (base_state, make_alias, BaseState::p0_comp , 1); MultiFab th0(base_state, make_alias, BaseState::th0_comp, 1); - if (anelastic == 1) { #ifdef _OPENMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) #endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); + for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + Box tbz = mfi.tilebox(); - // We don't compute a source term for z-momentum on the bottom or top boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); + // We don't compute a source term for z-momentum on the bottom or top boundary + if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); + if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); + const Array4 & cell_data = S_data[IntVars::cons].array(mfi); + const Array4 & cell_prim = S_prim.array(mfi); + const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - // Base state density and pressure - const Array4& r0_arr = r0.const_array(mfi); - const Array4& th0_arr = th0.const_array(mfi); + // Base state density and pressure + const Array4& r0_arr = r0.const_array(mfi); + const Array4& p0_arr = p0.const_array(mfi); + const Array4& th0_arr = th0.const_array(mfi); + + if ( anelastic && (solverChoice.moisture_type == MoistureType::None) ) + { + // ****************************************************************************************** + // Dry anelastic + // ****************************************************************************************** + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // + // Return -rho0 g (thetaprime / theta0) + // + buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k,grav_gpu[2], + r0_arr,th0_arr,cell_data); + }); + } + else if ( anelastic && (solverChoice.moisture_type != MoistureType::None) ) + { + // ****************************************************************************************** + // Moist anelastic + // ****************************************************************************************** + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + // + // Return -rho0 g (thetaprime / theta0) + // + buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k,grav_gpu[2],rv_over_rd, + r0_arr,th0_arr,cell_data); + }); + } + else if ( !anelastic && (solverChoice.moisture_type == MoistureType::None) ) + { + // ****************************************************************************************** + // Dry compressible + // ****************************************************************************************** + int n_q_dry = 0; + if (solverChoice.buoyancy_type == 1) { - if (solverChoice.moisture_type == MoistureType::None) { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // // Return -rho0 g (thetaprime / theta0) // - buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k, - grav_gpu[2], - r0_arr,th0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2], + r0_arr,cell_data); }); - } else { - // NOTE: For decomposition in the vertical direction, klo may not - // reside in the valid box and this call will yield an out - // of bounds error since it depends upon the surface theta_l + } + else if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 3) + { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // - // Return -rho0 g (thetaprime / theta0) + // Return -rho0 g (Tprime / T0) // - buoyancy_fab(i, j, k) = buoyancy_moist_anelastic(i,j,k, - grav_gpu[2],rv_over_rd, - r0_arr,th0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_dry_Tpert(i,j,k,grav_gpu[2],rd_over_cp, + r0_arr,p0_arr,th0_arr,cell_data); }); } - } // mfi - } - else - { - // ****************************************************************************************** - // Dry versions of buoyancy expressions (type 1 and type 2/3 -- types 2 and 3 are equivalent) - // ****************************************************************************************** - if (solverChoice.moisture_type == MoistureType::None) - { - int n_q_dry = 0; - if (solverChoice.buoyancy_type == 1) { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) - { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - // Base state density - const Array4& r0_arr = r0.const_array(mfi); - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data); - }); - } // mfi - - } - else // (buoyancy_type != 1) + else if (solverChoice.buoyancy_type == 4) { - // We now use the base state rather than planar average because - // 1) we don't want to average over the limited region of the fine level if doing multilevel. - // 2) it's cheaper to use the base state than to compute the horizontal averages - // 3) when running in a smallish domain, the horizontal average may evolve over time, - // which is not necessarily the intended behavior - // -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - // Base state density and pressure - const Array4& r0_arr = r0.const_array(mfi); - const Array4& p0_arr = p0.const_array(mfi); - const Array4& th0_arr = th0.const_array(mfi); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_dry_default(i,j,k, - grav_gpu[2],rd_over_cp, - r0_arr,p0_arr,th0_arr,cell_data); - }); - } // mfi - } // buoyancy_type - } // moisture type - else + // + // Return -rho0 g (Theta_prime / Theta_0) + // + buoyancy_fab(i, j, k) = buoyancy_dry_Thpert(i,j,k,grav_gpu[2], + r0_arr,th0_arr,cell_data); + }); + } // buoyancy_type for dry compressible + } + else // if ( !anelastic && (solverChoice.moisture_type != MoistureType::None) ) { - // ****************************************************************************************** - // Moist versions of buoyancy expressions - // ****************************************************************************************** + // ****************************************************************************************** + // Moist compressible + // ****************************************************************************************** - if ( (solverChoice.moisture_type == MoistureType::Kessler_NoRain) || - (solverChoice.moisture_type == MoistureType::SAM) || - (solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) ) - { - AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); - } - - if (solverChoice.buoyancy_type == 1) { - -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + if ( (solverChoice.moisture_type == MoistureType::Kessler_NoRain) || + (solverChoice.moisture_type == MoistureType::SAM) || + (solverChoice.moisture_type == MoistureType::SAM_NoPrecip_NoIce) ) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - // Base state density - const Array4& r0_arr = r0.const_array(mfi); + AMREX_ALWAYS_ASSERT(solverChoice.buoyancy_type == 1); + } + if (solverChoice.buoyancy_type == 1) + { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate, - grav_gpu[2],r0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_qstate,grav_gpu[2], + r0_arr,cell_data); }); - } // mfi - - } else { - - PlaneAverage state_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane); - PlaneAverage prim_ave(&S_prim , geom, solverChoice.ave_plane); - - // Compute horizontal averages of all components of each field - state_ave.compute_averages(ZDir(), state_ave.field()); - prim_ave.compute_averages(ZDir(), prim_ave.field()); - - int ncell = state_ave.ncell_line(); - - Gpu::HostVector rho_h(ncell), theta_h(ncell); - Gpu::DeviceVector rho_d(ncell), theta_d(ncell); - - state_ave.line_average(Rho_comp, rho_h); - Gpu::copyAsync(Gpu::hostToDevice, rho_h.begin(), rho_h.end(), rho_d.begin()); - - prim_ave.line_average(PrimTheta_comp, theta_h); - Gpu::copyAsync(Gpu::hostToDevice, theta_h.begin(), theta_h.end(), theta_d.begin()); - - Real* rho_d_ptr = rho_d.data(); - Real* theta_d_ptr = theta_d.data(); - - // Average valid moisture vars - Gpu::HostVector qv_h(ncell) , qc_h(ncell) , qp_h(ncell); - Gpu::DeviceVector qv_d(ncell,0.0), qc_d(ncell,0.0), qp_d(ncell,0.0); - if (n_qstate >=1) { - prim_ave.line_average(PrimQ1_comp, qv_h); - Gpu::copyAsync(Gpu::hostToDevice, qv_h.begin(), qv_h.end(), qv_d.begin()); - } - if (n_qstate >=2) { - prim_ave.line_average(PrimQ2_comp, qc_h); - Gpu::copyAsync(Gpu::hostToDevice, qc_h.begin(), qc_h.end(), qc_d.begin()); - } - if (n_qstate >=3) { - prim_ave.line_average(PrimQ3_comp, qp_h); - Gpu::copyAsync(Gpu::hostToDevice, qp_h.begin(), qp_h.end(), qp_d.begin()); } - Real* qv_d_ptr = qv_d.data(); - Real* qc_d_ptr = qc_d.data(); - Real* qp_d_ptr = qp_d.data(); - - if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 4 ) { + else if (solverChoice.buoyancy_type == 2 || solverChoice.buoyancy_type == 3) + { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4 & cell_prim = S_prim.array(mfi); - - // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp) - if (solverChoice.buoyancy_type == 2) { - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_type2(i,j,k,n_qstate,grav_gpu[2], - rho_d_ptr,theta_d_ptr, - qv_d_ptr,qc_d_ptr,qp_d_ptr, - cell_prim,cell_data); - }); - } else { - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_type4(i,j,k,n_qstate,grav_gpu[2], - rho_d_ptr,theta_d_ptr, - qv_d_ptr,qc_d_ptr,qp_d_ptr, - cell_prim,cell_data); - }); - } - } // mfi - - } else if (solverChoice.buoyancy_type == 3) { -#ifdef _OPENMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - for ( MFIter mfi(buoyancy,TilingIfNotGPU()); mfi.isValid(); ++mfi) + buoyancy_fab(i, j, k) = buoyancy_moist_Tpert(i,j,k,n_qstate,grav_gpu[2],rd_over_cp, + r0_arr,th0_arr,p0_arr, + cell_prim,cell_data); + }); + } + else if (solverChoice.buoyancy_type == 4) + { + ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - Box tbz = mfi.tilebox(); - - // We don't compute a source term for z-momentum on the bottom or top domain boundary - if (tbz.smallEnd(2) == klo) tbz.growLo(2,-1); - if (tbz.bigEnd(2) == khi) tbz.growHi(2,-1); - - const Array4< Real> & buoyancy_fab = buoyancy.array(mfi); - - const Array4 & cell_data = S_data[IntVars::cons].array(mfi); - const Array4 & cell_prim = S_prim.array(mfi); - - // TODO: ice has not been dealt with (q1=qv, q2=qv, q3=qp) - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - buoyancy_fab(i, j, k) = buoyancy_type3(i,j,k,n_qstate,grav_gpu[2], - rho_d_ptr,theta_d_ptr,qv_d_ptr, - cell_prim,cell_data); + buoyancy_fab(i, j, k) = buoyancy_moist_Thpert(i,j,k,n_qstate,grav_gpu[2], + r0_arr,th0_arr,cell_prim); }); - } // mfi - } // buoyancy_type - } // not buoyancy_type == 1 - } // has moisture - } // anelastic? + } + } // moist compressible + } // mfi } diff --git a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H index 499120896..00951a969 100644 --- a/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H +++ b/Source/TimeIntegration/ERF_TI_slow_rhs_fun.H @@ -15,6 +15,8 @@ // // Define primitive variables for all later RK stages // (We have already done this for the first RK step) + // Note that it is essential this happen before the call to make_mom_sources + // because some of the buoyancy routines use the primitive variables // if (nrk > 0) { int ng_cons = S_data[IntVars::cons].nGrow(); diff --git a/Submodules/AMReX b/Submodules/AMReX index f7f3ee9f5..1c8b54530 160000 --- a/Submodules/AMReX +++ b/Submodules/AMReX @@ -1 +1 @@ -Subproject commit f7f3ee9f57170e227af1fee1c6b95e5b41a45af8 +Subproject commit 1c8b5453023acb66bf65580b56e2f2f1dfd40c8e