From bc278b3b0afacf916722fe5459e7346eb78b5f07 Mon Sep 17 00:00:00 2001 From: "Aaron M. Lattanzi" <103702284+AMLattanzi@users.noreply.github.com> Date: Tue, 9 Jul 2024 12:04:18 -0700 Subject: [PATCH] BOMEX Subsidence mod (#1674) * Match FE source code. * Use FE qsat and balance with moisture. * Move wbar to z-faces and average. * Correct z location for wbar on z-face. * Revert qsat. --- Exec/DevTests/Bomex/prob.cpp | 8 +- Source/ERF.cpp | 2 +- Source/SourceTerms/ERF_make_mom_sources.cpp | 90 ++++++++++++--- Source/SourceTerms/ERF_make_sources.cpp | 118 ++++++++++++-------- Source/Utils/Microphysics_Utils.H | 36 ++++++ 5 files changed, 185 insertions(+), 69 deletions(-) diff --git a/Exec/DevTests/Bomex/prob.cpp b/Exec/DevTests/Bomex/prob.cpp index e6f2d286f..f368c7acf 100644 --- a/Exec/DevTests/Bomex/prob.cpp +++ b/Exec/DevTests/Bomex/prob.cpp @@ -325,7 +325,7 @@ Problem::update_w_subsidence (const Real& /*time*/, { if (wbar.empty()) return; - const int khi = geom.Domain().bigEnd()[2]; + const int khi = geom.Domain().bigEnd()[2] + 1; // lives on z-faces const Real* prob_lo = geom.ProbLo(); const auto dx = geom.CellSize(); @@ -343,9 +343,9 @@ Problem::update_w_subsidence (const Real& /*time*/, Real z_0 = 0.0; // (z_phys_cc) ? zlevels[0] : prob_lo[2] + 0.5 * dx[2]; Real slope1 = parms.wbar_sub_max / (parms.wbar_cutoff_max - z_0); Real slope2 = -parms.wbar_sub_max / (parms.wbar_cutoff_min - parms.wbar_cutoff_max); - // wbar[0] = 0.0; - for (int k = 0; k <= khi; k++) { - const Real z_cc = (z_phys_cc) ? zlevels[k] : prob_lo[2] + (k+0.5)* dx[2]; + wbar[0] = 0.0; + for (int k = 1; k <= khi; k++) { + const Real z_cc = (z_phys_cc) ? zlevels[k] : prob_lo[2] + k*dx[2]; if (z_cc <= parms.wbar_cutoff_max) { wbar[k] = slope1 * (z_cc - z_0); } else if (z_cc <= parms.wbar_cutoff_min) { diff --git a/Source/ERF.cpp b/Source/ERF.cpp index d0ba7d42f..2c00edfc1 100644 --- a/Source/ERF.cpp +++ b/Source/ERF.cpp @@ -774,7 +774,7 @@ ERF::InitData () h_w_subsid.resize(max_level+1, Vector(0)); d_w_subsid.resize(max_level+1, Gpu::DeviceVector(0)); for (int lev = 0; lev <= finest_level; lev++) { - const int domlen = geom[lev].Domain().length(2); + const int domlen = geom[lev].Domain().length(2) + 1; // lives on z-faces h_w_subsid[lev].resize(domlen, 0.0_rt); d_w_subsid[lev].resize(domlen, 0.0_rt); prob->update_w_subsidence(t_new[0], diff --git a/Source/SourceTerms/ERF_make_mom_sources.cpp b/Source/SourceTerms/ERF_make_mom_sources.cpp index 9dff15ee9..c7d72a2d5 100644 --- a/Source/SourceTerms/ERF_make_mom_sources.cpp +++ b/Source/SourceTerms/ERF_make_mom_sources.cpp @@ -110,12 +110,38 @@ void make_mom_sources (int /*level*/, // ***************************************************************************** // Planar averages for subsidence terms // ***************************************************************************** - Table1D dptr_u_plane, dptr_v_plane; - TableData u_plane_tab, v_plane_tab; + Table1D dptr_r_plane, dptr_u_plane, dptr_v_plane; + TableData r_plane_tab, u_plane_tab, v_plane_tab; if (dptr_wbar_sub) { - PlaneAverage u_ave(&xvel , geom, solverChoice.ave_plane, true); - PlaneAverage v_ave(&yvel , geom, solverChoice.ave_plane, true); + // Rho + PlaneAverage r_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true); + r_ave.compute_averages(ZDir(), r_ave.field()); + + int ncell = r_ave.ncell_line(); + Gpu::HostVector< Real> r_plane_h(ncell); + Gpu::DeviceVector< Real> r_plane_d(ncell); + + r_ave.line_average(Rho_comp, r_plane_h); + + Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin()); + + Real* dptr_r = r_plane_d.data(); + + IntVect ng_c = S_data[IntVars::cons].nGrowVect(); + Box tdomain = domain; tdomain.grow(2,ng_c[2]); + r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)}); + + int offset = ng_c[2]; + dptr_r_plane = r_plane_tab.table(); + ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept + { + dptr_r_plane(k-offset) = dptr_r[k]; + }); + + // U and V momentum + PlaneAverage u_ave(&(S_data[IntVars::xmom]), geom, solverChoice.ave_plane, true); + PlaneAverage v_ave(&(S_data[IntVars::ymom]), geom, solverChoice.ave_plane, true); u_ave.compute_averages(ZDir(), u_ave.field()); v_ave.compute_averages(ZDir(), v_ave.field()); @@ -125,8 +151,8 @@ void make_mom_sources (int /*level*/, Gpu::HostVector< Real> u_plane_h(u_ncell), v_plane_h(v_ncell); Gpu::DeviceVector< Real> u_plane_d(u_ncell), v_plane_d(v_ncell); - u_ave.line_average(0 , u_plane_h); - v_ave.line_average(0 , v_plane_h); + u_ave.line_average(0 , u_plane_h); + v_ave.line_average(0 , v_plane_h); Gpu::copyAsync(Gpu::hostToDevice, u_plane_h.begin(), u_plane_h.end(), u_plane_d.begin()); Gpu::copyAsync(Gpu::hostToDevice, v_plane_h.begin(), v_plane_h.end(), v_plane_d.begin()); @@ -134,8 +160,8 @@ void make_mom_sources (int /*level*/, Real* dptr_u = u_plane_d.data(); Real* dptr_v = v_plane_d.data(); - IntVect ng_u = xvel.nGrowVect(); - IntVect ng_v = yvel.nGrowVect(); + IntVect ng_u = S_data[IntVars::xmom].nGrowVect(); + IntVect ng_v = S_data[IntVars::ymom].nGrowVect(); Box udomain = domain; udomain.grow(2,ng_u[2]); Box vdomain = domain; vdomain.grow(2,ng_v[2]); u_plane_tab.resize({udomain.smallEnd(2)}, {udomain.bigEnd(2)}); @@ -280,16 +306,44 @@ void make_mom_sources (int /*level*/, // Add custom SUBSIDENCE terms // ***************************************************************************** if (solverChoice.custom_w_subsidence) { - ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - xmom_src_arr(i, j, k) -= dptr_wbar_sub[k] * - 0.5 * (dptr_u_plane(k+1) - dptr_u_plane(k-1)) * dxInv[2]; - }); - ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - ymom_src_arr(i, j, k) -= dptr_wbar_sub[k] * - 0.5 * (dptr_v_plane(k+1) - dptr_v_plane(k-1)) * dxInv[2]; - }); + if (solverChoice.custom_forcing_prim_vars) { + const int nr = Rho_comp; + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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]; + }); + ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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]; + }); + } else { + ParallelFor(tbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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]; + }); + ParallelFor(tby, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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]; + }); + } } // ***************************************************************************** diff --git a/Source/SourceTerms/ERF_make_sources.cpp b/Source/SourceTerms/ERF_make_sources.cpp index 69fbc5002..5b54a86bf 100644 --- a/Source/SourceTerms/ERF_make_sources.cpp +++ b/Source/SourceTerms/ERF_make_sources.cpp @@ -70,28 +70,50 @@ void make_sources (int level, // ***************************************************************************** // Planar averages for subsidence terms // ***************************************************************************** - Table1D dptr_t_plane, dptr_qv_plane, dptr_qc_plane; - TableData t_plane_tab, qv_plane_tab, qc_plane_tab; + Table1D dptr_r_plane, dptr_t_plane, dptr_qv_plane, dptr_qc_plane; + TableData r_plane_tab, t_plane_tab, qv_plane_tab, qc_plane_tab; if (dptr_wbar_sub) { - PlaneAverage t_ave(&S_prim, geom, solverChoice.ave_plane, true); + // Rho + PlaneAverage r_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true); + r_ave.compute_averages(ZDir(), r_ave.field()); + + int ncell = r_ave.ncell_line(); + Gpu::HostVector< Real> r_plane_h(ncell); + Gpu::DeviceVector< Real> r_plane_d(ncell); + + r_ave.line_average(Rho_comp, r_plane_h); + + Gpu::copyAsync(Gpu::hostToDevice, r_plane_h.begin(), r_plane_h.end(), r_plane_d.begin()); + + Real* dptr_r = r_plane_d.data(); + + IntVect ng_c = S_data[IntVars::cons].nGrowVect(); + Box tdomain = domain; tdomain.grow(2,ng_c[2]); + r_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)}); + + int offset = ng_c[2]; + dptr_r_plane = r_plane_tab.table(); + ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept + { + dptr_r_plane(k-offset) = dptr_r[k]; + }); + + // Rho * Theta + PlaneAverage t_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true); t_ave.compute_averages(ZDir(), t_ave.field()); - int ncell = t_ave.ncell_line(); Gpu::HostVector< Real> t_plane_h(ncell); Gpu::DeviceVector< Real> t_plane_d(ncell); - t_ave.line_average(PrimTheta_comp, t_plane_h); + t_ave.line_average(RhoTheta_comp, t_plane_h); Gpu::copyAsync(Gpu::hostToDevice, t_plane_h.begin(), t_plane_h.end(), t_plane_d.begin()); Real* dptr_t = t_plane_d.data(); - IntVect ng_c = S_prim.nGrowVect(); - Box tdomain = domain; tdomain.grow(2,ng_c[2]); t_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)}); - int offset = ng_c[2]; dptr_t_plane = t_plane_tab.table(); ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept { @@ -101,13 +123,13 @@ void make_sources (int level, if (solverChoice.moisture_type != MoistureType::None) { // Water vapor - PlaneAverage qv_ave(&S_prim, geom, solverChoice.ave_plane, true); + PlaneAverage qv_ave(&(S_data[IntVars::cons]), geom, solverChoice.ave_plane, true); qv_ave.compute_averages(ZDir(), qv_ave.field()); Gpu::HostVector< Real> qv_plane_h(ncell); Gpu::DeviceVector qv_plane_d(ncell); - qv_ave.line_average(PrimQ1_comp, qv_plane_h); + qv_ave.line_average(RhoQ1_comp, qv_plane_h); Gpu::copyAsync(Gpu::hostToDevice, qv_plane_h.begin(), qv_plane_h.end(), qv_plane_d.begin()); Real* dptr_qv = qv_plane_d.data(); @@ -119,26 +141,6 @@ void make_sources (int level, { dptr_qv_plane(k-offset) = dptr_qv[k]; }); - - // Cloud water - PlaneAverage qc_ave(&S_prim, geom, solverChoice.ave_plane, true); - qc_ave.compute_averages(ZDir(), qc_ave.field()); - - Gpu::HostVector< Real> qc_plane_h(ncell); - Gpu::DeviceVector qc_plane_d(ncell); - - qc_ave.line_average(PrimQ2_comp, qc_plane_h); - Gpu::copyAsync(Gpu::hostToDevice, qc_plane_h.begin(), qc_plane_h.end(), qc_plane_d.begin()); - - Real* dptr_qc = qc_plane_d.data(); - - qc_plane_tab.resize({tdomain.smallEnd(2)}, {tdomain.bigEnd(2)}); - - dptr_qc_plane = qc_plane_tab.table(); - ParallelFor(ncell, [=] AMREX_GPU_DEVICE (int k) noexcept - { - dptr_qc_plane(k-offset) = dptr_qc[k]; - }); } } @@ -239,12 +241,26 @@ void make_sources (int level, // ************************************************************************************* if (solverChoice.custom_w_subsidence) { const int n = RhoTheta_comp; - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - cell_src(i, j, k, n) -= dptr_wbar_sub[k] * - 0.5 * (dptr_t_plane(k+1) - dptr_t_plane(k-1)) * dxInv[2]; - - }); + if (solverChoice.custom_forcing_prim_vars) { + const int nr = Rho_comp; + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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]; + }); + } else { + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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]; + }); + } } // ************************************************************************************* @@ -253,16 +269,26 @@ void make_sources (int level, if (solverChoice.custom_w_subsidence && (solverChoice.moisture_type != MoistureType::None)) { const int nv = RhoQ1_comp; const int nc = RhoQ2_comp; - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - cell_src(i,j,k,nv) -= dptr_wbar_sub[k] * - 0.5 * (dptr_qv_plane(k+1) - dptr_qv_plane(k-1)) * dxInv[2]; - }); - ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - cell_src(i,j,k,nc) -= dptr_wbar_sub[k] * - 0.5 * (dptr_qc_plane(k+1) - dptr_qc_plane(k-1)) * dxInv[2]; - }); + if (solverChoice.custom_forcing_prim_vars) { + const int nr = Rho_comp; + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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 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]; + }); + } else { + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + 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 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]; + }); + } } // ************************************************************************************* diff --git a/Source/Utils/Microphysics_Utils.H b/Source/Utils/Microphysics_Utils.H index c5eb5b66f..1318221cd 100644 --- a/Source/Utils/Microphysics_Utils.H +++ b/Source/Utils/Microphysics_Utils.H @@ -42,6 +42,7 @@ amrex::Real erf_esati (amrex::Real t) { AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real erf_esatw (amrex::Real t) { +#if 1 amrex::Real const a0 = 6.105851; amrex::Real const a1 = 0.4440316; amrex::Real const a2 = 0.1430341e-1; @@ -62,6 +63,37 @@ amrex::Real erf_esatw (amrex::Real t) { esatw = 2.0*0.01*std::exp(9.550426 - 5723.265/t + 3.53068*std::log(t) - 0.00728332*t); } return esatw; +#else + amrex::Real a0 = 6.11239921; + amrex::Real a1 = 0.443987641; + amrex::Real a2 = 0.142986287e-1; + amrex::Real a3 = 0.264847430e-3; + amrex::Real a4 = 0.302950461e-5; + amrex::Real a5 = 0.206739458e-7; + amrex::Real a6 = 0.640689451e-10; + amrex::Real a7 = -0.952447341e-13; + amrex::Real a8 = -0.976195544e-15; + + amrex::Real a0i = 6.11147274; + amrex::Real a1i = 0.503160820; + amrex::Real a2i = 0.188439774e-1; + amrex::Real a3i = 0.420895665e-3; + amrex::Real a4i = 0.615021634e-5; + amrex::Real a5i = 0.602588177e-7; + amrex::Real a6i = 0.385852041e-9; + amrex::Real a7i = 0.146898966e-11; + amrex::Real a8i = 0.252751365e-14; + + amrex::Real Tc = std::max(-80.0, t-273.16); + + amrex::Real esatw; + if (t>=273.15){ + esatw = (a0 + Tc*(a1 +Tc*(a2 +Tc*(a3 +Tc*(a4 +Tc*(a5 +Tc*(a6 +Tc*(a7 +a8 *Tc))))))))*100.; + }else{ + esatw = (a0i + Tc*(a1i+Tc*(a2i+Tc*(a3i+Tc*(a4i+Tc*(a5i+Tc*(a6i+Tc*(a7i+a8i*Tc))))))))*100.; + } + return esatw; +#endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -122,7 +154,11 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void erf_qsatw (amrex::Real t, amrex::Real p, amrex::Real &qsatw) { amrex::Real esatw; esatw = erf_esatw(t); +#if 1 qsatw = Rd_on_Rv*esatw/std::max(esatw,p-esatw); +#else + qsatw = esatw/(R_v*t); +#endif } AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE