From cc6436e1913e7560ddf89bf3ca4c86481cb23252 Mon Sep 17 00:00:00 2001 From: AMLattanzi Date: Wed, 6 Nov 2024 09:17:07 -0800 Subject: [PATCH 1/2] Ran with kessler and sam successfully. --- Source/Microphysics/Kessler/ERF_Kessler.H | 3 + Source/Microphysics/Kessler/ERF_Kessler.cpp | 185 +++++++++++++------- Source/Microphysics/SAM/ERF_IceFall.cpp | 115 ++++++++---- Source/Microphysics/SAM/ERF_PrecipFall.cpp | 181 +++++++++++++------ Source/Microphysics/SAM/ERF_SAM.H | 3 + 5 files changed, 336 insertions(+), 151 deletions(-) diff --git a/Source/Microphysics/Kessler/ERF_Kessler.H b/Source/Microphysics/Kessler/ERF_Kessler.H index 1db7756fd..e1033a160 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.H +++ b/Source/Microphysics/Kessler/ERF_Kessler.H @@ -143,6 +143,9 @@ private: // Number of qstate variables int m_qstate_size = 3; + // CFL MAX for vertical advection + static constexpr amrex::Real CFL_MAX = 0.5; + // MicVar map (Qmoist indices -> MicVar enum) amrex::Vector MicVarMap; diff --git a/Source/Microphysics/Kessler/ERF_Kessler.cpp b/Source/Microphysics/Kessler/ERF_Kessler.cpp index 9906da3f8..5268f369a 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.cpp +++ b/Source/Microphysics/Kessler/ERF_Kessler.cpp @@ -22,53 +22,10 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) auto dm = tabs->DistributionMap(); fz.define(convert(ba, IntVect(0,0,1)), dm, 1, 0); // No ghost cells - Real dtn = dt; - - for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ - auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); - auto rain_accum_array = mic_fab_vars[MicVar_Kess::rain_accum]->array(mfi); - - auto fz_array = fz.array(mfi); - const Box& tbz = mfi.tilebox(); - - ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - Real rho_avg, qp_avg; - - if (k==k_lo) { - rho_avg = rho_array(i,j,k); - qp_avg = qp_array(i,j,k); - } else if (k==k_hi+1) { - rho_avg = rho_array(i,j,k-1); - qp_avg = qp_array(i,j,k-1); - } else { - rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3 - qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k)); - } - - qp_avg = std::max(0.0, qp_avg); - - Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s - - // NOTE: Fz is the sedimentation flux from the advective operator. - // In the terrain-following coordinate system, the z-deriv in - // the divergence uses the normal velocity (Omega). However, - // there are no u/v components to the sedimentation velocity. - // Therefore, we simply end up with a division by detJ when - // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv. - fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; - - if(k==k_lo){ - rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*qp_avg*V_terminal*dtn/1000.0*1000.0; // Divide by rho_water and convert to mm - } - - /*if(k==0){ - fz_array(i,j,k) = 0; - }*/ - }); - } + Real dtn = dt; + Real coef = dtn/dz; + // Saturation and evaporation calculations go first for ( MFIter mfi(*tabs,TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto qv_array = mic_fab_vars[MicVar_Kess::qv]->array(mfi); auto qc_array = mic_fab_vars[MicVar_Kess::qcl]->array(mfi); @@ -79,20 +36,13 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) auto theta_array = mic_fab_vars[MicVar_Kess::theta]->array(mfi); auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); - const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4{}; - - const auto& box3d = mfi.tilebox(); - - auto fz_array = fz.array(mfi); + const auto& tbx = mfi.tilebox(); // Expose for GPU Real d_fac_cond = m_fac_cond; - ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - // Jacobian determinant - Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0; - qv_array(i,j,k) = std::max(0.0, qv_array(i,j,k)); qc_array(i,j,k) = std::max(0.0, qc_array(i,j,k)); qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); @@ -163,14 +113,9 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) dq_clwater_to_rain = std::min(dq_clwater_to_rain, qc_array(i,j,k)); } - if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; - if(std::fabs(fz_array(i,j,k )) < 1e-14) fz_array(i,j,k ) = 0.0; - Real dq_sed = dtn * dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k))/dz; - if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0; - qv_array(i,j,k) += -dq_vapor_to_clwater + dq_clwater_to_vapor + dq_rain_to_vapor; qc_array(i,j,k) += dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain; - qp_array(i,j,k) += dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; + qp_array(i,j,k) += dq_clwater_to_rain - dq_rain_to_vapor; Real theta_over_T = theta_array(i,j,k)/tabs_array(i,j,k); theta_array(i,j,k) += theta_over_T * d_fac_cond * (dq_vapor_to_clwater - dq_clwater_to_vapor - dq_rain_to_vapor); @@ -182,8 +127,126 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) qt_array(i,j,k) = qv_array(i,j,k) + qc_array(i,j,k); }); } + + // Precompute terminal velocity for substepping + for ( MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi ){ + auto rho_array = mic_fab_vars[MicVar_Kess::rho]->array(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + + auto fz_array = fz.array(mfi); + const Box& tbz = mfi.tilebox(); + + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real rho_avg, qp_avg; + + if (k==k_lo) { + rho_avg = rho_array(i,j,k); + qp_avg = qp_array(i,j,k); + } else if (k==k_hi+1) { + rho_avg = rho_array(i,j,k-1); + qp_avg = qp_array(i,j,k-1); + } else { + rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3 + qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k)); + } + + qp_avg = std::max(0.0, qp_avg); + + Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s + + // NOTE: Fz is the sedimentation flux from the advective operator. + // In the terrain-following coordinate system, the z-deriv in + // the divergence uses the normal velocity (Omega). However, + // there are no u/v components to the sedimentation velocity. + // Therefore, we simply end up with a division by detJ when + // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv. + fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; + }); + } // mfi + + // Compute number of substeps from maximum terminal velocity + Real wt_max; + int n_substep; + auto const& ma_fz_arr = fz.const_arrays(); + GpuTuple max = ParReduce(TypeList{}, + TypeList{}, + fz, IntVect(0), + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept + -> GpuTuple + { + return { ma_fz_arr[box_no](i,j,k) }; + }); + wt_max = get<0>(max) + 1.0e-16; + n_substep = int( std::ceil(wt_max * coef / CFL_MAX) ); + AMREX_ALWAYS_ASSERT(n_substep >= 1); + coef /= Real(n_substep); + dtn /= Real(n_substep); + + // Substep the vertical advection + for (int nsub(0); nsubarray(mfi); + auto qp_array = mic_fab_vars[MicVar_Kess::qp]->array(mfi); + auto rain_accum_array = mic_fab_vars[MicVar_Kess::rain_accum]->array(mfi); + auto fz_array = fz.array(mfi); + + const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4{}; + + const Box& tbx = mfi.tilebox(); + const Box& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0)); + + // Update vertical flux every substep + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real rho_avg, qp_avg; + if (k==k_lo) { + rho_avg = rho_array(i,j,k); + qp_avg = qp_array(i,j,k); + } else if (k==k_hi+1) { + rho_avg = rho_array(i,j,k-1); + qp_avg = qp_array(i,j,k-1); + } else { + rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); // Convert to g/cm^3 + qp_avg = 0.5*(qp_array(i,j,k-1) + qp_array(i,j,k)); + } + + qp_avg = std::max(0.0, qp_avg); + + Real V_terminal = 36.34*std::pow(rho_avg*0.001*qp_avg, 0.1346)*std::pow(rho_avg/1.16, -0.5); // in m/s + + // NOTE: Fz is the sedimentation flux from the advective operator. + // In the terrain-following coordinate system, the z-deriv in + // the divergence uses the normal velocity (Omega). However, + // there are no u/v components to the sedimentation velocity. + // Therefore, we simply end up with a division by detJ when + // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv. + fz_array(i,j,k) = rho_avg*V_terminal*qp_avg; + + if(k==k_lo){ + rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*qp_avg*V_terminal*dtn/1000.0*1000.0; // Divide by rho_water and convert to mm + } + }); + + // Update precip every substep + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Jacobian determinant + Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0; + + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; + if(std::fabs(fz_array(i,j,k )) < 1e-14) fz_array(i,j,k ) = 0.0; + Real dq_sed = dJinv * (1.0/rho_array(i,j,k)) * (fz_array(i,j,k+1) - fz_array(i,j,k)) * coef; + if(std::fabs(dq_sed) < 1e-14) dq_sed = 0.0; + + qp_array(i,j,k) += dq_sed; + qp_array(i,j,k) = std::max(0.0, qp_array(i,j,k)); + }); + } // mfi + } // nsub } + if (solverChoice.moisture_type == MoistureType::Kessler_NoRain){ // get the temperature, dentisy, theta, qt and qc from input diff --git a/Source/Microphysics/SAM/ERF_IceFall.cpp b/Source/Microphysics/SAM/ERF_IceFall.cpp index 9518a9401..31f80a195 100644 --- a/Source/Microphysics/SAM/ERF_IceFall.cpp +++ b/Source/Microphysics/SAM/ERF_IceFall.cpp @@ -1,6 +1,7 @@ #include #include "ERF_SAM.H" #include "ERF_TileNoZ.H" +#include using namespace amrex; @@ -55,7 +56,7 @@ void SAM::IceFall (const SolverChoice& sc) { rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); qci_avg = 0.5*(qci_array(i,j,k-1) + qci_array(i,j,k)); } - Real vt_ice = min( 0.4 , 8.66 * pow( (max(0.,qci_avg)+1.e-10) , 0.24) ); + Real vt_ice = min( 0.4 , 8.66 * pow( (std::max(0.,qci_avg)+1.e-10) , 0.24) ); // NOTE: Fz is the sedimentation flux from the advective operator. // In the terrain-following coordinate system, the z-deriv in @@ -67,37 +68,85 @@ void SAM::IceFall (const SolverChoice& sc) { }); } - for (MFIter mfi(*qci, TileNoZ()); mfi.isValid(); ++mfi) { - auto qci_array = qci->array(mfi); - auto qn_array = qn->array(mfi); - auto qt_array = qt->array(mfi); - auto rho_array = rho->array(mfi); - auto fz_array = fz.array(mfi); - - const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4{}; - - const auto& box3d = mfi.tilebox(); - - ParallelFor(box3d, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Jacobian determinant - Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0; - - //================================================== - // Cloud ice sedimentation (A32) - //================================================== - Real dqi = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef; - dqi = std::max(-qci_array(i,j,k), dqi); - - // Add this increment to both non-precipitating and total water. - qci_array(i,j,k) += dqi; - qn_array(i,j,k) += dqi; - qt_array(i,j,k) += dqi; - - // NOTE: Sedimentation does not affect the potential temperature, - // but it does affect the liquid/ice static energy. - // No source to Theta occurs here. - }); - } + // Compute number of substeps from maximum terminal velocity + Real wt_max; + int n_substep; + auto const& ma_fz_arr = fz.const_arrays(); + GpuTuple max = ParReduce(TypeList{}, + TypeList{}, + fz, IntVect(0), + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept + -> GpuTuple + { + return { ma_fz_arr[box_no](i,j,k) }; + }); + wt_max = get<0>(max) + 1.0e-16; + n_substep = int( std::ceil(wt_max * coef / CFL_MAX) ); + AMREX_ALWAYS_ASSERT(n_substep >= 1); + coef /= Real(n_substep); + dtn /= Real(n_substep); + + // Substep the vertical advection + for (int nsub(0); nsubarray(mfi); + auto qn_array = qn->array(mfi); + auto qt_array = qt->array(mfi); + auto rho_array = rho->array(mfi); + auto fz_array = fz.array(mfi); + + const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4{}; + + const auto& tbx = mfi.tilebox(); + const auto& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0)); + + // Update vertical flux every substep + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real rho_avg, qci_avg; + if (k==k_lo) { + rho_avg = rho_array(i,j,k); + qci_avg = qci_array(i,j,k); + } else if (k==k_hi+1) { + rho_avg = rho_array(i,j,k-1); + qci_avg = qci_array(i,j,k-1); + } else { + rho_avg = 0.5*(rho_array(i,j,k-1) + rho_array(i,j,k)); + qci_avg = 0.5*(qci_array(i,j,k-1) + qci_array(i,j,k)); + } + Real vt_ice = min( 0.4 , 8.66 * pow( (std::max(0.,qci_avg)+1.e-10) , 0.24) ); + + // NOTE: Fz is the sedimentation flux from the advective operator. + // In the terrain-following coordinate system, the z-deriv in + // the divergence uses the normal velocity (Omega). However, + // there are no u/v components to the sedimentation velocity. + // Therefore, we simply end up with a division by detJ when + // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv. + fz_array(i,j,k) = rho_avg*vt_ice*qci_avg; + }); + + // Update precip every substep + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Jacobian determinant + Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0; + + //================================================== + // Cloud ice sedimentation (A32) + //================================================== + Real dqi = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef; + dqi = std::max(-qci_array(i,j,k), dqi); + + // Add this increment to both non-precipitating and total water. + qci_array(i,j,k) += dqi; + qn_array(i,j,k) += dqi; + qt_array(i,j,k) += dqi; + + // NOTE: Sedimentation does not affect the potential temperature, + // but it does affect the liquid/ice static energy. + // No source to Theta occurs here. + }); + } // mfi + } // nsub } diff --git a/Source/Microphysics/SAM/ERF_PrecipFall.cpp b/Source/Microphysics/SAM/ERF_PrecipFall.cpp index abbfec8fd..0590d8990 100644 --- a/Source/Microphysics/SAM/ERF_PrecipFall.cpp +++ b/Source/Microphysics/SAM/ERF_PrecipFall.cpp @@ -1,6 +1,7 @@ #include "ERF_Constants.H" #include "ERF_SAM.H" #include "ERF_TileNoZ.H" +#include using namespace amrex; @@ -57,15 +58,12 @@ SAM::PrecipFall (const SolverChoice& sc) SAM_moisture_type = 2; } - // Add sedimentation of precipitation field to the vert. vel. + // Precompute the vertical fluxes for CFL constraint for (MFIter mfi(fz, TilingIfNotGPU()); mfi.isValid(); ++mfi) { auto qp_array = qp->array(mfi); auto rho_array = rho->array(mfi); auto tabs_array = tabs->array(mfi); auto fz_array = fz.array(mfi); - auto rain_accum_array = rain_accum->array(mfi); - auto snow_accum_array = snow_accum->array(mfi); - auto graup_accum_array = graup_accum->array(mfi); const auto& box3d = mfi.tilebox(); @@ -111,66 +109,135 @@ SAM::PrecipFall (const SolverChoice& sc) // Therefore, we simply end up with a division by detJ when // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv. fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg); + }); + } + + // Compute number of substeps from maximum terminal velocity + Real wt_max; + int n_substep; + auto const& ma_fz_arr = fz.const_arrays(); + GpuTuple max = ParReduce(TypeList{}, + TypeList{}, + fz, IntVect(0), + [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k) noexcept + -> GpuTuple + { + return { ma_fz_arr[box_no](i,j,k) }; + }); + wt_max = get<0>(max) + 1.0e-16; + n_substep = int( std::ceil(wt_max * coef / CFL_MAX) ); + AMREX_ALWAYS_ASSERT(n_substep >= 1); + coef /= Real(n_substep); + dtn /= Real(n_substep); + + // Substep the vertical advection + for (int nsub(0); nsubarray(mfi); + auto qps_array = qps->array(mfi); + auto qpg_array = qpg->array(mfi); + auto qp_array = qp->array(mfi); + auto rho_array = rho->array(mfi); + auto tabs_array = tabs->array(mfi); + auto fz_array = fz.array(mfi); + + auto rain_accum_array = rain_accum->array(mfi); + auto snow_accum_array = snow_accum->array(mfi); + auto graup_accum_array = graup_accum->array(mfi); + + const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4{}; + + const auto& tbx = mfi.tilebox(); + const auto& tbz = mfi.tilebox(IntVect(0,0,1),IntVect(0)); + + // Update vertical flux every substep + ParallelFor(tbz, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + Real rho_avg, tab_avg, qp_avg; + if (k==k_lo) { + rho_avg = rho_array(i,j,k); + tab_avg = tabs_array(i,j,k); + qp_avg = qp_array(i,j,k); + } else if (k==k_hi+1) { + rho_avg = rho_array(i,j,k-1); + tab_avg = tabs_array(i,j,k-1); + qp_avg = qp_array(i,j,k-1); + } else { + rho_avg = 0.5*( rho_array(i,j,k-1) + rho_array(i,j,k)); + tab_avg = 0.5*(tabs_array(i,j,k-1) + tabs_array(i,j,k)); + qp_avg = 0.5*( qp_array(i,j,k-1) + qp_array(i,j,k)); + } + + Real Pprecip = 0.0; + if(qp_avg > qp_threshold) { + Real omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + } + Real qrr = omp*qp_avg; + Real qss = (1.0-omp)*(1.0-omg)*qp_avg; + Real qgg = (1.0-omp)*(omg)*qp_avg; + Pprecip = omp*vrain*std::pow(rho_avg*qrr,1.0+crain) + + (1.0-omp)*( (1.0-omg)*vsnow*std::pow(rho_avg*qss,1.0+csnow) + + omg *vgrau*std::pow(rho_avg*qgg,1.0+cgrau) ); + } - if(k==k_lo){ + // NOTE: Fz is the sedimentation flux from the advective operator. + // In the terrain-following coordinate system, the z-deriv in + // the divergence uses the normal velocity (Omega). However, + // there are no u/v components to the sedimentation velocity. + // Therefore, we simply end up with a division by detJ when + // evaluating the source term: dJinv * (flux_hi - flux_lo) * dzinv. + fz_array(i,j,k) = Pprecip * std::sqrt(rho_0/rho_avg); + + if(k==k_lo){ + Real omp, omg; + if (SAM_moisture_type == 2) { + omp = 1.0; + omg = 0.0; + } else { + omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + } + rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*(omp*qp_avg)*vrain*dtn/rhor*1000.0; // Divide by rho_water and convert to mm + snow_accum_array(i,j,k) = snow_accum_array(i,j,k) + rho_avg*(1.0-omp)*(1.0-omg)*qp_avg*vrain*dtn/rhos*1000.0; // Divide by rho_snow and convert to mm + graup_accum_array(i,j,k) = graup_accum_array(i,j,k) + rho_avg*(1.0-omp)*(omg)*qp_avg*vrain*dtn/rhog*1000.0; // Divide by rho_graupel and convert to mm + } + }); + + // Update precip every substep + ParallelFor(tbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Jacobian determinant + Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0; + + //================================================== + // Precipitating sedimentation (A19) + //================================================== + Real dqp = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef; Real omp, omg; if (SAM_moisture_type == 2) { omp = 1.0; omg = 0.0; } else { - omp = std::max(0.0,std::min(1.0,(tab_avg-tprmin)*a_pr)); - omg = std::max(0.0,std::min(1.0,(tab_avg-tgrmin)*a_gr)); + omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); + omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); } - rain_accum_array(i,j,k) = rain_accum_array(i,j,k) + rho_avg*(omp*qp_avg)*vrain*dtn/rhor*1000.0; // Divide by rho_water and convert to mm - snow_accum_array(i,j,k) = snow_accum_array(i,j,k) + rho_avg*(1.0-omp)*(1.0-omg)*qp_avg*vrain*dtn/rhos*1000.0; // Divide by rho_snow and convert to mm - graup_accum_array(i,j,k) = graup_accum_array(i,j,k) + rho_avg*(1.0-omp)*(omg)*qp_avg*vrain*dtn/rhog*1000.0; // Divide by rho_graupel and convert to mm - } - - }); - } - - for (MFIter mfi(*qp, TileNoZ()); mfi.isValid(); ++mfi) { - auto qpr_array = qpr->array(mfi); - auto qps_array = qps->array(mfi); - auto qpg_array = qpg->array(mfi); - auto qp_array = qp->array(mfi); - auto rho_array = rho->array(mfi); - auto tabs_array = tabs->array(mfi); - auto fz_array = fz.array(mfi); - - const auto dJ_array = (m_detJ_cc) ? m_detJ_cc->const_array(mfi) : Array4{}; - - const auto& box3d = mfi.tilebox(); - // Update precipitation mass fraction and liquid-ice static - // energy using precipitation fluxes computed in this column. - ParallelFor(box3d, [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - // Jacobian determinant - Real dJinv = (dJ_array) ? 1.0/dJ_array(i,j,k) : 1.0; - - //================================================== - // Precipitating sedimentation (A19) - //================================================== - Real dqp = dJinv * (1.0/rho_array(i,j,k)) * ( fz_array(i,j,k+1) - fz_array(i,j,k) ) * coef; - Real omp, omg; - if (SAM_moisture_type == 2) { - omp = 1.0; - omg = 0.0; - } else { - omp = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tprmin)*a_pr)); - omg = std::max(0.0,std::min(1.0,(tabs_array(i,j,k)-tgrmin)*a_gr)); - } - - qpr_array(i,j,k) = std::max(0.0, qpr_array(i,j,k) + dqp*omp); - qps_array(i,j,k) = std::max(0.0, qps_array(i,j,k) + dqp*(1.0-omp)*(1.0-omg)); - qpg_array(i,j,k) = std::max(0.0, qpg_array(i,j,k) + dqp*(1.0-omp)*omg); - qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); - - // NOTE: Sedimentation does not affect the potential temperature, - // but it does affect the liquid/ice static energy. - // No source to Theta occurs here. - }); - } // mfi + qpr_array(i,j,k) = std::max(0.0, qpr_array(i,j,k) + dqp*omp); + qps_array(i,j,k) = std::max(0.0, qps_array(i,j,k) + dqp*(1.0-omp)*(1.0-omg)); + qpg_array(i,j,k) = std::max(0.0, qpg_array(i,j,k) + dqp*(1.0-omp)*omg); + qp_array(i,j,k) = qpr_array(i,j,k) + qps_array(i,j,k) + qpg_array(i,j,k); + + // NOTE: Sedimentation does not affect the potential temperature, + // but it does affect the liquid/ice static energy. + // No source to Theta occurs here. + }); + } // mfi + } // nsub } diff --git a/Source/Microphysics/SAM/ERF_SAM.H b/Source/Microphysics/SAM/ERF_SAM.H index d2e1d0929..e137129e8 100644 --- a/Source/Microphysics/SAM/ERF_SAM.H +++ b/Source/Microphysics/SAM/ERF_SAM.H @@ -286,6 +286,9 @@ private: // Number of qmoist variables int m_qstate_size = 6; + // CFL MAX for vertical advection + static constexpr amrex::Real CFL_MAX = 0.5; + // MicVar map (Qmoist indices -> MicVar enum) amrex::Vector MicVarMap; From a0aa61502a6bd8cad72ef17e17ff96c66e6a22d6 Mon Sep 17 00:00:00 2001 From: Aaron Lattanzi Date: Wed, 6 Nov 2024 12:27:45 -0800 Subject: [PATCH 2/2] Use numerical limits. --- Source/Microphysics/Kessler/ERF_Kessler.cpp | 2 +- Source/Microphysics/SAM/ERF_IceFall.cpp | 2 +- Source/Microphysics/SAM/ERF_PrecipFall.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/Microphysics/Kessler/ERF_Kessler.cpp b/Source/Microphysics/Kessler/ERF_Kessler.cpp index 5268f369a..958eb6440 100644 --- a/Source/Microphysics/Kessler/ERF_Kessler.cpp +++ b/Source/Microphysics/Kessler/ERF_Kessler.cpp @@ -177,7 +177,7 @@ void Kessler::AdvanceKessler (const SolverChoice &solverChoice) { return { ma_fz_arr[box_no](i,j,k) }; }); - wt_max = get<0>(max) + 1.0e-16; + wt_max = get<0>(max) + std::numeric_limits::epsilon(); n_substep = int( std::ceil(wt_max * coef / CFL_MAX) ); AMREX_ALWAYS_ASSERT(n_substep >= 1); coef /= Real(n_substep); diff --git a/Source/Microphysics/SAM/ERF_IceFall.cpp b/Source/Microphysics/SAM/ERF_IceFall.cpp index 31f80a195..99bca7f52 100644 --- a/Source/Microphysics/SAM/ERF_IceFall.cpp +++ b/Source/Microphysics/SAM/ERF_IceFall.cpp @@ -80,7 +80,7 @@ void SAM::IceFall (const SolverChoice& sc) { { return { ma_fz_arr[box_no](i,j,k) }; }); - wt_max = get<0>(max) + 1.0e-16; + wt_max = get<0>(max) + std::numeric_limits::epsilon(); n_substep = int( std::ceil(wt_max * coef / CFL_MAX) ); AMREX_ALWAYS_ASSERT(n_substep >= 1); coef /= Real(n_substep); diff --git a/Source/Microphysics/SAM/ERF_PrecipFall.cpp b/Source/Microphysics/SAM/ERF_PrecipFall.cpp index 0590d8990..07daef1d8 100644 --- a/Source/Microphysics/SAM/ERF_PrecipFall.cpp +++ b/Source/Microphysics/SAM/ERF_PrecipFall.cpp @@ -124,7 +124,7 @@ SAM::PrecipFall (const SolverChoice& sc) { return { ma_fz_arr[box_no](i,j,k) }; }); - wt_max = get<0>(max) + 1.0e-16; + wt_max = get<0>(max) + std::numeric_limits::epsilon(); n_substep = int( std::ceil(wt_max * coef / CFL_MAX) ); AMREX_ALWAYS_ASSERT(n_substep >= 1); coef /= Real(n_substep);