From 23d6ef921a2f23eefef0f2545f4825f50565f681 Mon Sep 17 00:00:00 2001 From: Ann Almgren Date: Wed, 20 Nov 2024 15:28:14 -0800 Subject: [PATCH] Modify dry compressible buoyancy to use thetapert not Tpert --- Source/SourceTerms/ERF_buoyancy_utils.H | 155 +++++++++++++---------- Source/SourceTerms/ERF_make_buoyancy.cpp | 19 ++- 2 files changed, 97 insertions(+), 77 deletions(-) diff --git a/Source/SourceTerms/ERF_buoyancy_utils.H b/Source/SourceTerms/ERF_buoyancy_utils.H index e8db1015e..9f73da597 100644 --- a/Source/SourceTerms/ERF_buoyancy_utils.H +++ b/Source/SourceTerms/ERF_buoyancy_utils.H @@ -7,55 +7,34 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real -buoyancy_dry_anelastic (int& i, +buoyancy_dry_thetapert (int& i, int& j, int& k, amrex::Real const& grav_gpu, const amrex::Array4& r0_arr, - const amrex::Array4& p0_arr, + const amrex::Array4& th0_arr, const amrex::Array4& cell_data) { + // + // This returns - g rho0 (thetaprime / theta) + // // Note: this is the same term as the moist anelastic buoyancy when qv = qc = qt = 0 - amrex::Real theta_d_lo = cell_data(i,j,k-1,RhoTheta_comp)/r0_arr(i,j,k-1); - amrex::Real theta_d_hi = cell_data(i,j,k ,RhoTheta_comp)/r0_arr(i,j,k); - - amrex::Real theta_d_wface = amrex::Real(0.5) * (theta_d_lo + theta_d_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_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_0_wface = amrex::Real(0.5) * (theta_d_0_lo + theta_d_0_hi); - - return (-grav_gpu * (theta_d_wface - theta_d_0_wface) / theta_d_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1))); -} + // We don't do divide by 2 here because the "2" will cancel below + amrex::Real theta_d_zface = theta_d_lo + theta_d_hi; -AMREX_GPU_DEVICE -AMREX_FORCE_INLINE -amrex::Real -buoyancy_dry_anelastic_T (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& 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; + // We don't do divide by 2 here because the "2" will cancel below + amrex::Real theta_d_0_zface = th0_arr(i,j,k) + th0_arr(i,j,k-1); - 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; + // rho0 on z-face + amrex::Real r0_zface = 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 -grav_gpu * r0_zface * (theta_d_zface - theta_d_0_zface) / theta_d_0_zface; } - AMREX_GPU_DEVICE AMREX_FORCE_INLINE amrex::Real @@ -68,6 +47,9 @@ buoyancy_moist_anelastic (int& i, const amrex::Array4& th0_arr, const amrex::Array4& cell_data) { + // + // This returns - g rho0 (thetaprime / theta) + // 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); @@ -90,41 +72,17 @@ buoyancy_moist_anelastic (int& i, return (-grav_gpu * (theta_v_wface - theta_v_0_wface) / theta_v_0_wface * 0.5 * (r0_arr(i,j,k) + r0_arr(i,j,k-1))); } -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 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 -buoyancy_type1 (int& i, - int& j, - int& k, - const int& n_qstate, - amrex::Real const& grav_gpu, - const amrex::Array4& r0_arr, - const amrex::Array4& cell_data) +buoyancy_rhopert (int& i, + int& j, + int& k, + const int& n_qstate, + amrex::Real const& grav_gpu, + const amrex::Array4& r0_arr, + const amrex::Array4& 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); @@ -272,4 +230,67 @@ buoyancy_type4 (int& i, return (-qavg * r0avg * grav_gpu); } + +// ******************************************************************************** +// The functions below this are not currently in use +// ******************************************************************************** + +#if 0 +AMREX_GPU_DEVICE +AMREX_FORCE_INLINE +amrex::Real +buoyancy_dry_Tpert_comp (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) +{ + // + // This returns - g rho0 (Tprime / T) + // + 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 +buoyancy_dry_Tpert_anel (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& cell_data) +{ + // + // This returns - g rho0 (Tprime / T) + // + 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); +} +#endif #endif diff --git a/Source/SourceTerms/ERF_make_buoyancy.cpp b/Source/SourceTerms/ERF_make_buoyancy.cpp index ca101605e..6b8dfcd55 100644 --- a/Source/SourceTerms/ERF_make_buoyancy.cpp +++ b/Source/SourceTerms/ERF_make_buoyancy.cpp @@ -70,15 +70,14 @@ void make_buoyancy (Vector& S_data, // 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 (solverChoice.moisture_type == MoistureType::None) { ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - buoyancy_fab(i, j, k) = buoyancy_dry_anelastic(i,j,k, + buoyancy_fab(i, j, k) = buoyancy_dry_thetapert(i,j,k, grav_gpu[2], - r0_arr,p0_arr,cell_data); + r0_arr,th0_arr,cell_data); }); } else { // NOTE: For decomposition in the vertical direction, klo may not @@ -121,7 +120,7 @@ void make_buoyancy (Vector& S_data, ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - buoyancy_fab(i, j, k) = buoyancy_type1(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data); + buoyancy_fab(i, j, k) = buoyancy_rhopert(i,j,k,n_q_dry,grav_gpu[2],r0_arr,cell_data); }); } // mfi @@ -147,7 +146,7 @@ void make_buoyancy (Vector& S_data, // Base state density and pressure const Array4& r0_arr = r0.const_array(mfi); - const Array4& p0_arr = p0.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); @@ -155,9 +154,9 @@ void make_buoyancy (Vector& S_data, 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); + buoyancy_fab(i, j, k) = buoyancy_dry_thetapert(i,j,k, + grav_gpu[2], + r0_arr,th0_arr,cell_data); }); } // mfi } // buoyancy_type @@ -196,8 +195,8 @@ void make_buoyancy (Vector& S_data, ParallelFor(tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - buoyancy_fab(i, j, k) = buoyancy_type1(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