diff --git a/src/BulkMicrophysicsTendencies.jl b/src/BulkMicrophysicsTendencies.jl index f5ed437d7..be4b90d39 100644 --- a/src/BulkMicrophysicsTendencies.jl +++ b/src/BulkMicrophysicsTendencies.jl @@ -283,90 +283,92 @@ Returns a `NamedTuple` containing the nonzero entries of `M` and `e`. e2 = zero(FT) e4 = zero(FT) - # --- Phase change: vapor ↔ cloud condensate --- - D = src.S_phase_change_vap_lcl / max(q_min, q_lcl) - is_source = src.S_phase_change_vap_lcl >= zero(FT) - e1 += ifelse(is_source, src.S_phase_change_vap_lcl, zero(FT)) - M11 += ifelse(is_source, zero(FT), D) - - D = src.S_phase_change_vap_icl / max(q_min, q_icl) - is_source = src.S_phase_change_vap_icl >= zero(FT) - e2 += ifelse(is_source, src.S_phase_change_vap_icl, zero(FT)) - M22 += ifelse(is_source, zero(FT), D) - - # --- Melt: ice cloud → liquid cloud --- - D = src.S_melt_icl_lcl / max(q_min, q_icl) - M22 -= D - M12 += D - - # --- Autoconversion: donor-based transfer --- - D = src.S_acnv_lcl_rai / max(q_min, q_lcl) - M11 -= D - M31 += D - - D = src.S_acnv_icl_sno / max(q_min, q_icl) - M22 -= D - M42 += D - - # --- Accretion: donor-based transfer --- - D = src.S_accr_lcl_rai / max(q_min, q_lcl) - M11 -= D - M31 += D - - # lcl + sno accretion (cold/warm arms already zeroed) - D_cold = src.S_accr_lcl_sno_cold / max(q_min, q_lcl) - D_warm = src.S_accr_lcl_sno_warm / max(q_min, q_lcl) - M11 -= D_cold + D_warm - M31 += D_warm # warm: lcl → rai - M41 += D_cold # cold: lcl → sno - - # thermal melt of sno from warm lcl - D = src.S_accr_melt_lcl_sno / max(q_min, q_sno) - M44 -= D - M34 += D - - D = src.S_accr_icl_rai / max(q_min, q_icl) - M22 -= D - M42 += D - - D = src.S_accr_icl_sno / max(q_min, q_icl) - M22 -= D - M42 += D - - # rain frozen in icl + rai collision - D = src.S_accr_freeze_icl_rai / max(q_min, q_rai) - M33 -= D - M43 += D - - # warm arm: sno melts → rai (already zero when cold) - D = src.S_accr_rai_sno_warm / max(q_min, q_sno) - M44 -= D - M34 += D - - # thermal melt of sno from warm rai (already zero when cold) - D = src.S_accr_melt_rai_sno / max(q_min, q_sno) - M44 -= D - M34 += D - - # cold arm: rai freezes → sno (already zero when warm) - D = src.S_accr_rai_sno_cold / max(q_min, q_rai) - M33 -= D - M43 += D - - # --- Rain phase change: sink to vapor (always zero or negative) --- - D = (-src.S_phase_change_vap_rai) / max(q_min, q_rai) - M33 -= D - - # --- Snow phase change: deposition/sublimation --- - D = src.S_phase_change_vap_sno / max(q_min, q_sno) - is_source = src.S_phase_change_vap_sno >= zero(FT) - e4 += ifelse(is_source, src.S_phase_change_vap_sno, zero(FT)) - M44 += ifelse(is_source, zero(FT), D) - - # --- Snow melt: snow → rain --- - D = src.S_melt_sno_rai / max(q_min, q_sno) - M44 -= D - M34 += D + @fastmath begin + # --- Phase change: vapor ↔ cloud condensate --- + D = src.S_phase_change_vap_lcl / max(q_min, q_lcl) + is_source = src.S_phase_change_vap_lcl >= zero(FT) + e1 += ifelse(is_source, src.S_phase_change_vap_lcl, zero(FT)) + M11 += ifelse(is_source, zero(FT), D) + + D = src.S_phase_change_vap_icl / max(q_min, q_icl) + is_source = src.S_phase_change_vap_icl >= zero(FT) + e2 += ifelse(is_source, src.S_phase_change_vap_icl, zero(FT)) + M22 += ifelse(is_source, zero(FT), D) + + # --- Melt: ice cloud → liquid cloud --- + D = src.S_melt_icl_lcl / max(q_min, q_icl) + M22 -= D + M12 += D + + # --- Autoconversion: donor-based transfer --- + D = src.S_acnv_lcl_rai / max(q_min, q_lcl) + M11 -= D + M31 += D + + D = src.S_acnv_icl_sno / max(q_min, q_icl) + M22 -= D + M42 += D + + # --- Accretion: donor-based transfer --- + D = src.S_accr_lcl_rai / max(q_min, q_lcl) + M11 -= D + M31 += D + + # lcl + sno accretion (cold/warm arms already zeroed) + D_cold = src.S_accr_lcl_sno_cold / max(q_min, q_lcl) + D_warm = src.S_accr_lcl_sno_warm / max(q_min, q_lcl) + M11 -= D_cold + D_warm + M31 += D_warm # warm: lcl → rai + M41 += D_cold # cold: lcl → sno + + # thermal melt of sno from warm lcl + D = src.S_accr_melt_lcl_sno / max(q_min, q_sno) + M44 -= D + M34 += D + + D = src.S_accr_icl_rai / max(q_min, q_icl) + M22 -= D + M42 += D + + D = src.S_accr_icl_sno / max(q_min, q_icl) + M22 -= D + M42 += D + + # rain frozen in icl + rai collision + D = src.S_accr_freeze_icl_rai / max(q_min, q_rai) + M33 -= D + M43 += D + + # warm arm: sno melts → rai (already zero when cold) + D = src.S_accr_rai_sno_warm / max(q_min, q_sno) + M44 -= D + M34 += D + + # thermal melt of sno from warm rai (already zero when cold) + D = src.S_accr_melt_rai_sno / max(q_min, q_sno) + M44 -= D + M34 += D + + # cold arm: rai freezes → sno (already zero when warm) + D = src.S_accr_rai_sno_cold / max(q_min, q_rai) + M33 -= D + M43 += D + + # --- Rain phase change: sink to vapor (always zero or negative) --- + D = (-src.S_phase_change_vap_rai) / max(q_min, q_rai) + M33 -= D + + # --- Snow phase change: deposition/sublimation --- + D = src.S_phase_change_vap_sno / max(q_min, q_sno) + is_source = src.S_phase_change_vap_sno >= zero(FT) + e4 += ifelse(is_source, src.S_phase_change_vap_sno, zero(FT)) + M44 += ifelse(is_source, zero(FT), D) + + # --- Snow melt: snow → rain --- + D = src.S_melt_sno_rai / max(q_min, q_sno) + M44 -= D + M34 += D + end return ( M11 = M11, M12 = M12, M22 = M22, @@ -392,8 +394,14 @@ The system uses a sparse structure specific to the 1-moment microphysics model. Because sinks are linearized as `-D q`, they are effectively integrated as exponential decays over the substep. + +`@noinline` is a deliberate GPU register-pressure device-function barrier: this function +returns only 4 Float32, which keeps the call cheap, while preventing the ClimaAtmos.jl SGS +quadrature evaluator (which calls this ~9 times per grid column) from inlining all 9 copies +into one giant kernel. Measured -26% device time (255→154 regs) on the reference longrun +config — do not remove without re-profiling. """ -@inline function _linearized_implicit_step( +@noinline function _linearized_implicit_step( ::Microphysics1Moment, mp::CMP.Microphysics1MParams, tps, ρ, T, q_tot, q_lcl, q_icl, q_rai, q_sno, Δt, ) @@ -407,46 +415,48 @@ exponential decays over the substep. q_min = TDI.TD.Parameters.q_min(tps) lin = _linearize(src, q_lcl, q_icl, q_rai, q_sno, q_min) - invΔt = one(FT) / Δt - - # A = I/Δt - M - a11 = invΔt - lin.M11 - a12 = -lin.M12 - a22 = invΔt - lin.M22 - a31 = -lin.M31 - a33 = invΔt - lin.M33 - a34 = -lin.M34 - a41 = -lin.M41 - a42 = -lin.M42 - a43 = -lin.M43 - a44 = invΔt - lin.M44 - - # rhs = e + q_0/Δt - # e3 = 0 by the 1m model - b1 = lin.e1 + invΔt * q_lcl - b2 = lin.e2 + invΔt * q_icl - b3 = invΔt * q_rai - b4 = lin.e4 + invΔt * q_sno - - # Solve 2×2 system for q_lcl, q_icl (coupled via ice melt M12) - det12 = a11 * a22 # a21 = 0 - q_lcl_new = (b1 * a22 - a12 * b2) / det12 - q_icl_new = a11 * b2 / det12 - - # Reduced 2x2 system for q_rai_new, q_sno_new - r3 = muladd(-a31, q_lcl_new, b3) - r4 = muladd(-a41, q_lcl_new, muladd(-a42, q_icl_new, b4)) - - det = muladd(-a34, a43, a33 * a44) - # det is a positive number because a44 and a33 are positive (greater than invΔt) - # and a34 and a43 are non-positive so we don't need to safeguard division by det. - q_rai_new = (r3 * a44 - a34 * r4) / det - q_sno_new = (a33 * r4 - r3 * a43) / det - - dq_lcl_dt = (q_lcl_new - q_lcl) * invΔt - dq_icl_dt = (q_icl_new - q_icl) * invΔt - dq_rai_dt = (q_rai_new - q_rai) * invΔt - dq_sno_dt = (q_sno_new - q_sno) * invΔt + @fastmath begin + invΔt = one(FT) / Δt + + # A = I/Δt - M + a11 = invΔt - lin.M11 + a12 = -lin.M12 + a22 = invΔt - lin.M22 + a31 = -lin.M31 + a33 = invΔt - lin.M33 + a34 = -lin.M34 + a41 = -lin.M41 + a42 = -lin.M42 + a43 = -lin.M43 + a44 = invΔt - lin.M44 + + # rhs = e + q_0/Δt + # e3 = 0 by the 1m model + b1 = lin.e1 + invΔt * q_lcl + b2 = lin.e2 + invΔt * q_icl + b3 = invΔt * q_rai + b4 = lin.e4 + invΔt * q_sno + + # Solve 2×2 system for q_lcl, q_icl (coupled via ice melt M12) + det12 = a11 * a22 # a21 = 0 + q_lcl_new = (b1 * a22 - a12 * b2) / det12 + q_icl_new = a11 * b2 / det12 + + # Reduced 2x2 system for q_rai_new, q_sno_new + r3 = muladd(-a31, q_lcl_new, b3) + r4 = muladd(-a41, q_lcl_new, muladd(-a42, q_icl_new, b4)) + + det = muladd(-a34, a43, a33 * a44) + # det is a positive number because a44 and a33 are positive (greater than invΔt) + # and a34 and a43 are non-positive so we don't need to safeguard division by det. + q_rai_new = (r3 * a44 - a34 * r4) / det + q_sno_new = (a33 * r4 - r3 * a43) / det + + dq_lcl_dt = (q_lcl_new - q_lcl) * invΔt + dq_icl_dt = (q_icl_new - q_icl) * invΔt + dq_rai_dt = (q_rai_new - q_rai) * invΔt + dq_sno_dt = (q_sno_new - q_sno) * invΔt + end return (; dq_lcl_dt, dq_icl_dt, dq_rai_dt, dq_sno_dt) end diff --git a/src/Common.jl b/src/Common.jl index 87877ea28..aa1c5c28a 100644 --- a/src/Common.jl +++ b/src/Common.jl @@ -58,8 +58,8 @@ numerical robustness. D_vapor_safe = max(D_vapor, UT.ϵ_numerics(FT)) K_therm_safe = max(K_therm, UT.ϵ_numerics(FT)) - return 1 / - (L / K_therm_safe / T * (L / R_v / T - 1) + R_v * T / D_vapor_safe / p_vs_safe) + return @fastmath 1 / + (L / K_therm_safe / T * (L / R_v / T - 1) + R_v * T / D_vapor_safe / p_vs_safe) end """ @@ -97,8 +97,8 @@ numerical robustness. D_vapor_safe = max(D_vapor, UT.ϵ_numerics(FT)) K_therm_safe = max(K_therm, UT.ϵ_numerics(FT)) - return 1 / - (L / K_therm_safe / T * (L / R_v / T - 1) + R_v * T / D_vapor_safe / p_vs_safe) + return @fastmath 1 / + (L / K_therm_safe / T * (L / R_v / T - 1) + R_v * T / D_vapor_safe / p_vs_safe) end """ @@ -163,10 +163,10 @@ This curve smoothly transitions from y = 0 for 0 < x < x_0 to y = x - x_0 for x_ # translation of the curve in x and y to enforce zero at x = 0 # Using log1mexp for numerical stability: log1mexp(x) = log(1 - exp(x)) - trnslt = -LEF.log1mexp(-k) / k - kt = k * (x_safe / x_0_safe - 1 + trnslt) + trnslt = @fastmath -LEF.log1mexp(-k) / k + kt = @fastmath k * (x_safe / x_0_safe - 1 + trnslt) # log1pexp handles all kt values correctly - result = (LEF.log1pexp(kt) / k - trnslt) * x_0_safe + result = @fastmath (LEF.log1pexp(kt) / k - trnslt) * x_0_safe # Handle edge cases: if x ≈ 0, return 0; if x_0 ≈ 0, return x return ifelse(x < UT.ϵ_numerics(FT), FT(0), ifelse(x_0 < UT.ϵ_numerics(FT), x, result)) diff --git a/src/Microphysics1M.jl b/src/Microphysics1M.jl index b0bfdc646..d26cf4bb8 100644 --- a/src/Microphysics1M.jl +++ b/src/Microphysics1M.jl @@ -82,7 +82,7 @@ Returns the intercept parameter of the assumed Marshall-Palmer distribution """ @inline function get_n0((; ν, μ)::CMP.ParticlePDFSnow{FT}, q_sno::FT, ρ::FT) where {FT} safe_q_sno = max(q_sno, UT.ϵ_numerics(FT)) - return ifelse(q_sno > UT.ϵ_numerics(FT), μ * (ρ * safe_q_sno)^ν, zero(FT)) + return ifelse(q_sno > UT.ϵ_numerics(FT), (@fastmath μ * (ρ * safe_q_sno)^ν), zero(FT)) end @inline get_n0((; n0)::CMP.ParticlePDFIceRain{FT}, args...) where {FT} = n0 @@ -99,9 +99,11 @@ Guards against unphysical density ratios (ρ > ρw) that would cause sqrt of neg - `ρ`: air density (rain only) """ @inline function get_v0((; C_drag, ρw, grav, r0)::CMP.Blk1MVelTypeRain{FT}, ρ::FT) where {FT} - # Guard against ρ > ρw (unphysical but could occur from numerical errors) - density_factor = max(ρw / ρ - 1, zero(FT)) - return sqrt(FT(8 / 3) / C_drag * density_factor * grav * r0) + @fastmath begin + # Guard against ρ > ρw (unphysical but could occur from numerical errors) + density_factor = max(ρw / ρ - 1, zero(FT)) + return sqrt(FT(8 / 3) / C_drag * density_factor * grav * r0) + end end @inline get_v0((; v0)::CMP.Blk1MVelTypeSnow{FT}, args...) where {FT} = v0 @@ -146,8 +148,10 @@ average particles. The value is clipped at `r0 * 1e-5` to prevent numerical issu # ParticleMass constructor for GPU performance. qp = UT.clamp_to_nonneg(q) ρp = UT.clamp_to_nonneg(ρ) - denom = χm * m0 * max(n0, UT.ϵ_numerics(FT)) * gamma_coeff - λ_inv = (ρp * qp * r0^(me + Δm) / denom)^(1 / (me + Δm + 1)) + @fastmath begin + denom = χm * m0 * max(n0, UT.ϵ_numerics(FT)) * gamma_coeff + λ_inv = (ρp * qp * r0^(me + Δm) / denom)^(1 / (me + Δm + 1)) + end return max(r0 * FT(1e-5), λ_inv) end @@ -233,7 +237,7 @@ Fall velocity of individual particles is parameterized: # gamma_term = SF.gamma(me + ve + Δm + Δv + 1) (pre-computed in vel) # gamma_coeff = SF.gamma(me + Δm + 1) (pre-computed in mass) - fall_w = χv * v0 * (λ_inv / r0)^(ve + Δv) * gamma_term / gamma_coeff + fall_w = @fastmath χv * v0 * (λ_inv / r0)^(ve + Δv) * gamma_term / gamma_coeff return ifelse(q > UT.ϵ_numerics(FT), fall_w, zero(FT)) end @@ -354,7 +358,7 @@ using the prescribed cloud droplet number concentration. @inline function conv_q_lcl_to_q_rai(opt::CMP.Kessler1M, mp, tps, micro, thermo) q_lcl = micro.q_lcl (; τ, q_threshold, k) = opt.acnv1M - return CO.logistic_function_integral(q_lcl, q_threshold, k) / τ + return @fastmath CO.logistic_function_integral(q_lcl, q_threshold, k) / τ end @inline function conv_q_lcl_to_q_rai(opt::CMP.PrescribedNd, mp, tps, micro, thermo) @@ -414,7 +418,7 @@ Harrington et al. (1995) and Kaul et al. (2015). @inline function conv_q_icl_to_q_sno(opt::CMP.NoSupersaturation, mp, tps, micro, thermo, sd = nothing) (; τ, q_threshold, k) = opt.acnv1M q_icl = micro.q_icl - return CO.logistic_function_integral(q_icl, q_threshold, k) / τ + return @fastmath CO.logistic_function_integral(q_icl, q_threshold, k) / τ end @inline function conv_q_icl_to_q_sno( @@ -461,7 +465,7 @@ contribution of warm liquid on snow. T_freeze = TDI.T_freeze(tps) ΔT = T - T_freeze is_cold = (T <= T_freeze) - return ifelse(is_cold, zero(T), cv_l / L_f * ΔT) + return ifelse(is_cold, zero(T), @fastmath(cv_l / L_f * ΔT)) end # Gating convention for the process-rate kernels below: the rate is computed @@ -506,8 +510,8 @@ Internal low-level kernel. Prefer the option-dispatched API. # gamma_accr = SF.gamma(ae + ve + Δa + Δv + 1) (pre-computed in vel) accr_rate = - q_clo * E * n0 * a0 * v0 * χa * χv * λ_inv * - gamma_accr / (r0 / λ_inv)^(ae + ve + Δa + Δv) + @fastmath q_clo * E * n0 * a0 * v0 * χa * χv * λ_inv * + gamma_accr / (r0 / λ_inv)^(ae + ve + Δa + Δv) cond = q_clo > UT.ϵ_numerics(FT) && q_pre > UT.ϵ_numerics(FT) return ifelse(cond, accr_rate, zero(FT)) @@ -552,9 +556,9 @@ end # gamma_accr_rain_sink = SF.gamma(me + ae + ve + Δm + Δa + Δv + 1) (pre-computed in vel) accr_rate = - E / ρ * n0 * n0_ice * m0 * a0 * v0 * χm * χa * χv * λ_ice_inv * λ_inv * - gamma_accr_rain_sink / - (r0 / λ_inv)^FT(me + ae + ve + Δm + Δa + Δv) + @fastmath E / ρ * n0 * n0_ice * m0 * a0 * v0 * χm * χa * χv * λ_ice_inv * λ_inv * + gamma_accr_rain_sink / + (r0 / λ_inv)^FT(me + ae + ve + Δm + Δa + Δv) cond = q_icl > UT.ϵ_numerics(FT) && q_rai > UT.ϵ_numerics(FT) return ifelse(cond, accr_rate, zero(FT)) @@ -624,20 +628,22 @@ deviations are proportional to the mean fall velocities, with coefficient v_ti = terminal_velocity(type_i, blk1mveltype_ti, ρ, q_i, v0_i, λ_i_inv) v_tj = terminal_velocity(type_j, blk1mveltype_tj, ρ, q_j, v0_j, λ_j_inv) - # Add simple parameterization for velocity dispersion, assuming that fall velocity - # standard deviations are proportional to the mean fall velocities, with coefficient - # coeff_disp - Δv_eff = sqrt((v_ti - v_tj)^2 + coeff_disp * (v_ti^2 + v_tj^2)) - - # We use the recurrence relation Γ(x+1) = xΓ(x) to simplify gamma terms. - # gamma_coeff = Γ(δ + 1) is pre-computed. - accr_rate = - FT(π) / ρ * n0_i * n0_j * m0 * χm * E_ij * Δv_eff * gamma_coeff / - r0^δ * ( - 2 * λ_i_inv^3 * λ_j_inv^(δ + 1) + - 2 * (δ + 1) * λ_i_inv^2 * λ_j_inv^(δ + 2) + - (δ + 2) * (δ + 1) * λ_i_inv * λ_j_inv^(δ + 3) - ) + @fastmath begin + # Add simple parameterization for velocity dispersion, assuming that fall velocity + # standard deviations are proportional to the mean fall velocities, with coefficient + # coeff_disp + Δv_eff = sqrt((v_ti - v_tj)^2 + coeff_disp * (v_ti^2 + v_tj^2)) + + # We use the recurrence relation Γ(x+1) = xΓ(x) to simplify gamma terms. + # gamma_coeff = Γ(δ + 1) is pre-computed. + accr_rate = + FT(π) / ρ * n0_i * n0_j * m0 * χm * E_ij * Δv_eff * gamma_coeff / + r0^δ * ( + 2 * λ_i_inv^3 * λ_j_inv^(δ + 1) + + 2 * (δ + 1) * λ_i_inv^2 * λ_j_inv^(δ + 2) + + (δ + 2) * (δ + 1) * λ_i_inv * λ_j_inv^(δ + 3) + ) + end cond = q_i > UT.ϵ_numerics(FT) && q_j > UT.ϵ_numerics(FT) return ifelse(cond, accr_rate, zero(FT)) @@ -942,17 +948,19 @@ Only evaporation is considered (sub-saturated over liquid); result is clamped a_vent = vent.a b_vent = vent.b - Sc = ν_air / max(D_vapor, UT.ϵ_numerics(FT)) + @fastmath begin + Sc = ν_air / max(D_vapor, UT.ϵ_numerics(FT)) - evap_rate = - 4 * FT(π) * n0 / ρ * S * G * λ_inv^2 * - ( - a_vent + - b_vent * cbrt(Sc) / - (r0 / λ_inv)^((ve + Δv) / 2) * - sqrt(2 * v0 * χv / ν_air * λ_inv) * - gamma_vent - ) + evap_rate = + 4 * FT(π) * n0 / ρ * S * G * λ_inv^2 * + ( + a_vent + + b_vent * cbrt(Sc) / + (r0 / λ_inv)^((ve + Δv) / 2) * + sqrt(2 * v0 * χv / ν_air * λ_inv) * + gamma_vent + ) + end cond = q_rai > UT.ϵ_numerics(FT) && S < FT(0) return min(zero(FT), ifelse(cond, evap_rate, zero(FT))) @@ -1019,17 +1027,19 @@ end a_vent = vent.a b_vent = vent.b - Sc = ν_air / max(D_vapor, UT.ϵ_numerics(FT)) + @fastmath begin + Sc = ν_air / max(D_vapor, UT.ϵ_numerics(FT)) - subl_rate = - 4 * FT(π) * n0 / ρ * S * G * λ_inv^2 * - ( - a_vent + - b_vent * cbrt(Sc) / - (r0 / λ_inv)^((ve + Δv) / 2) * - sqrt(2 * v0 * χv / ν_air * λ_inv) * - gamma_vent - ) + subl_rate = + 4 * FT(π) * n0 / ρ * S * G * λ_inv^2 * + ( + a_vent + + b_vent * cbrt(Sc) / + (r0 / λ_inv)^((ve + Δv) / 2) * + sqrt(2 * v0 * χv / ν_air * λ_inv) * + gamma_vent + ) + end cond = q_sno > UT.ϵ_numerics(FT) return ifelse(cond, subl_rate, zero(FT)) @@ -1069,7 +1079,7 @@ Returns the tendency due to cloud ice melt. L = TDI.Lf(tps, T) (; n0) = pdf λ_inv = sd.λ_inv_icl - cloud_ice_melt_rate = 4 * FT(π) * n0 / ρ * K_therm / L * (T - T_freeze) * λ_inv^2 + cloud_ice_melt_rate = @fastmath 4 * FT(π) * n0 / ρ * K_therm / L * (T - T_freeze) * λ_inv^2 cond = q_icl > UT.ϵ_numerics(FT) && T > T_freeze return ifelse(cond, cloud_ice_melt_rate, zero(FT)) @@ -1120,18 +1130,20 @@ Returns the tendency due to snow melt. a_vent = vent.a b_vent = vent.b - # Schmidt number (guard against division by near-zero D_vapor) - Sc = ν_air / max(D_vapor, UT.ϵ_numerics(FT)) - - snow_melt_rate = - 4 * FT(π) * n0 / ρ * K_therm / L * (T - T_freeze) * λ_inv^2 * - ( - a_vent + - b_vent * cbrt(Sc) / - (r0 / λ_inv)^((ve + Δv) / 2) * - sqrt(2 * v0 * χv / ν_air * λ_inv) * - gamma_vent - ) + @fastmath begin + # Schmidt number (guard against division by near-zero D_vapor) + Sc = ν_air / max(D_vapor, UT.ϵ_numerics(FT)) + + snow_melt_rate = + 4 * FT(π) * n0 / ρ * K_therm / L * (T - T_freeze) * λ_inv^2 * + ( + a_vent + + b_vent * cbrt(Sc) / + (r0 / λ_inv)^((ve + Δv) / 2) * + sqrt(2 * v0 * χv / ν_air * λ_inv) * + gamma_vent + ) + end cond = q_sno > UT.ϵ_numerics(FT) && T > T_freeze return ifelse(cond, snow_melt_rate, zero(FT)) diff --git a/src/MicrophysicsNonEq.jl b/src/MicrophysicsNonEq.jl index e4619f0c5..027763677 100644 --- a/src/MicrophysicsNonEq.jl +++ b/src/MicrophysicsNonEq.jl @@ -72,7 +72,7 @@ temperature for a given phase of water. - `T` - temperature [K] """ @inline function dqcld_dT(qᵥ_sat, L, Rᵥ, T) - return qᵥ_sat * (L / (Rᵥ * T^2) - 1 / T) + return @fastmath qᵥ_sat * (L / (Rᵥ * T^2) - 1 / T) end """ @@ -86,7 +86,7 @@ Computes the thermodynamic adjustment factor Γ. - `dqcld_dT` - derivative of saturation specific humidity with respect to temperature [kg/kg/K] """ @inline function gamma_helper(L, cₚ_air, dqcld_dT) - return 1 + (L / cₚ_air) * dqcld_dT + return @fastmath 1 + (L / cₚ_air) * dqcld_dT end """ @@ -128,7 +128,7 @@ Morrison & Milbrandt (2015), https://doi.org/10.1175/JAS-D-14-0065.1. timescale = τ * Γₗ # compute the tendency - return ifelse( + return @fastmath ifelse( sat_excess < 0, -min(-sat_excess, max(0, q_lcl)) / timescale, sat_excess / timescale, @@ -175,7 +175,7 @@ Morrison & Milbrandt (2015), https://doi.org/10.1175/JAS-D-14-0065.1. timescale = τ * Γᵢ # compute the tendency - tendency = ifelse( + tendency = @fastmath ifelse( sat_excess < 0, -min(-sat_excess, max(0, q_icl)) / timescale, sat_excess / timescale, diff --git a/test/ad_compat_tests.jl b/test/ad_compat_tests.jl index 03909c18b..d433612aa 100644 --- a/test/ad_compat_tests.jl +++ b/test/ad_compat_tests.jl @@ -131,9 +131,9 @@ function test_ad_compatibility(FT) J_fd = similar(J) for j in 1:8 h = 1e-6 * r.x[j] - xp = copy(r.x); + xp = copy(r.x) xp[j] += h - xm = copy(r.x); + xm = copy(r.x) xm[j] -= h J_fd[:, j] = (f(xp) - f(xm)) / 2h end diff --git a/test/gpu_performance.jl b/test/gpu_performance.jl index 44a3c0782..ce0e7a904 100644 --- a/test/gpu_performance.jl +++ b/test/gpu_performance.jl @@ -159,17 +159,17 @@ function run_gpu_performance_benchmarks(FT) end end - a = FT(2.5); - x = FT(3.0); - p = FT(0.6); + a = FT(2.5) + x = FT(3.0) + p = FT(0.6) q = FT(0.4) b_sf = BT.@benchmark SF.gamma_inc($a, $x) b_ut = BT.@benchmark UT.gamma_inc($a, $x) - @info "gamma_inc benchmark (SF vs UT):" SF=BT.minimum(b_sf) UT=BT.minimum(b_ut) + @info "gamma_inc benchmark (SF vs UT):" SF = BT.minimum(b_sf) UT = BT.minimum(b_ut) b_inv_sf = BT.@benchmark SF.gamma_inc_inv($a, $p, $q) b_inv_ut = BT.@benchmark UT.gamma_inc_inv($a, $p, $q) - @info "gamma_inc_inv benchmark (SF vs UT):" SF=BT.minimum(b_inv_sf) UT=BT.minimum(b_inv_ut) + @info "gamma_inc_inv benchmark (SF vs UT):" SF = BT.minimum(b_inv_sf) UT = BT.minimum(b_inv_ut) end # 2. Kernel Benchmarks (1-Moment BMT, 2-Moment BMT, P3) diff --git a/test/p3_rho_d_stability.jl b/test/p3_rho_d_stability.jl index 00d2388a3..84671cf32 100644 --- a/test/p3_rho_d_stability.jl +++ b/test/p3_rho_d_stability.jl @@ -28,6 +28,7 @@ end @test ρ_g ≈ Float32(ρ_g_ref(BigFloat(F_rim), BigFloat(ρ_rim))) rtol = 1.0f-5 end # The value and its derivative with respect to F_rim stay finite under ForwardDiff in Float32. - g(x) = (ρ_d = P3.get_ρ_d(mass, x, 4.0f2); P3.get_ρ_g(x, 4.0f2, ρ_d)) + g(x) = (ρ_d = P3.get_ρ_d(mass, x, 4.0f2); + P3.get_ρ_g(x, 4.0f2, ρ_d)) @test isfinite(FD.derivative(g, 1.0f-4)) end