diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index b7fecefa54..52f6adc196 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -198,7 +198,10 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) advect_tke = use_prognostic_tke(turbconv_model) point_type = eltype(Fields.coordinate_field(Y.c)) (; dt) = p - ᶜJ = Fields.local_geometry_field(Y.c).J + ᶜJ, ᶠJ = Fields.local_geometry_field(Y.c).J, Fields.local_geometry_field(Y.f).J + ᶜρ, ᶠρ = Y.c.ρ, face_density(Y.c.ρ) + ᶠρJ = @. lazy(ᶠρ * ᶠJ) + ᶜρJ = @. lazy(ᶜρ * ᶜJ) (; ᶜf³, ᶠf¹², ᶜΦ) = p.core (; ᶜu, ᶠu³, ᶜK, ᶜts) = p.precomputed (; edmfx_mse_q_tot_upwinding) = n > 0 || advect_tke ? p.atmos.numerics : all_nothing @@ -217,14 +220,14 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) advect_tke ? ( turbconv_model isa PrognosticEDMFX ? - (@. lazy(ρa⁰(Y.c.ρ, Y.c.sgsʲs, turbconv_model))) : Y.c.ρ + (@. lazy(ρa⁰(ᶜρ, Y.c.sgsʲs, turbconv_model))) : ᶜρ ) : nothing ᶜρ⁰ = if advect_tke if n > 0 (; ᶜts⁰) = p.precomputed @. lazy(TD.air_density(thermo_params, ᶜts⁰)) else - Y.c.ρ + ᶜρ end else nothing @@ -254,14 +257,12 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) end # Without the CT12(), the right-hand side would be a CT1 or CT2 in 2D space. - ᶜρ = Y.c.ρ - # Full vertical advection of passive tracers (like liq, rai, etc) ... # If sgs_mass_flux is true, the advection term is computed from the sum of SGS fluxes if p.atmos.edmfx_model.sgs_mass_flux isa Val{false} foreach_gs_tracer(Yₜ, Y) do ᶜρχₜ, ᶜρχ, ρχ_name if !(ρχ_name in (@name(ρe_tot), @name(ρq_tot))) - ᶜχ = @. lazy(specific(ᶜρχ, Y.c.ρ)) + ᶜχ = @. lazy(specific(ᶜρχ, ᶜρ)) vtt = vertical_transport(ᶜρ, ᶠu³, ᶜχ, FT(dt), tracer_upwinding) @. ᶜρχₜ += vtt end @@ -283,7 +284,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) end if !(p.atmos.moisture_model isa DryModel) && energy_q_tot_upwinding != Val(:none) - ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) + ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, ᶜρ)) vtt = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, FT(dt), energy_q_tot_upwinding) vtt_central = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, FT(dt), Val(:none)) @. Yₜ.c.ρq_tot += vtt - vtt_central @@ -291,9 +292,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) if isnothing(ᶠf¹²) # shallow atmosphere - @. Yₜ.c.uₕ -= - ᶜinterp(ᶠω¹² × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / (Y.c.ρ * ᶜJ) + - (ᶜf³ + ᶜω³) × CT12(ᶜu) + @. Yₜ.c.uₕ -= ᶜinterp(ᶠω¹² × (ᶠρJ * ᶠu³)) / ᶜρJ + (ᶜf³ + ᶜω³) × CT12(ᶜu) @. Yₜ.f.u₃ -= ᶠω¹² × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK) for j in 1:n @. Yₜ.f.sgsʲs.:($$j).u₃ -= @@ -302,9 +301,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t) end else # deep atmosphere - @. Yₜ.c.uₕ -= - ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / - (Y.c.ρ * ᶜJ) + (ᶜf³ + ᶜω³) × CT12(ᶜu) + @. Yₜ.c.uₕ -= ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠρJ * ᶠu³)) / ᶜρJ + (ᶜf³ + ᶜω³) × CT12(ᶜu) @. Yₜ.f.u₃ -= (ᶠf¹² + ᶠω¹²) × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK) for j in 1:n @. Yₜ.f.sgsʲs.:($$j).u₃ -= @@ -349,11 +346,7 @@ Modifies EDMF updraft fields in `Yₜ.c.sgsʲs` and `Yₜ.f.sgsʲs`. """ edmfx_sgs_vertical_advection_tendency!(Yₜ, Y, p, t, turbconv_model) = nothing -function edmfx_sgs_vertical_advection_tendency!( - Yₜ, - Y, - p, - t, +function edmfx_sgs_vertical_advection_tendency!(Yₜ, Y, p, t, turbconv_model::PrognosticEDMFX, ) (; params) = p @@ -395,26 +388,21 @@ function edmfx_sgs_vertical_advection_tendency!( end for j in 1:n - ᶜa = (@. lazy(draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j)))) + ᶜρʲ = ᶜρʲs.:($j) + ᶜa = (@. lazy(draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲ))) # Flux form vertical advection of area farction with the grid mean velocity vtt = - vertical_transport(ᶜρʲs.:($j), ᶠu³ʲs.:($j), ᶜa, dt, edmfx_mse_q_tot_upwinding) + vertical_transport(ᶜρʲ, ᶠu³ʲs.:($j), ᶜa, dt, edmfx_mse_q_tot_upwinding) @. Yₜ.c.sgsʲs.:($$j).ρa += vtt # Advective form advection of mse and q_tot with the grid mean velocity # Note: This allocates because the function is too long - va = vertical_advection( - ᶠu³ʲs.:($j), - Y.c.sgsʲs.:($j).mse, - edmfx_mse_q_tot_upwinding, - ) + va = vertical_advection(ᶠu³ʲs.:($j), Y.c.sgsʲs.:($j).mse, edmfx_mse_q_tot_upwinding) @. Yₜ.c.sgsʲs.:($$j).mse += va va = vertical_advection( - ᶠu³ʲs.:($j), - Y.c.sgsʲs.:($j).q_tot, - edmfx_mse_q_tot_upwinding, + ᶠu³ʲs.:($j), Y.c.sgsʲs.:($j).q_tot, edmfx_mse_q_tot_upwinding ) @. Yₜ.c.sgsʲs.:($$j).q_tot += va @@ -433,15 +421,8 @@ function edmfx_sgs_vertical_advection_tendency!( ᶜ∂ρ∂t_sed = p.scratch.ᶜtemp_scalar_3 @. ᶜ∂ρ∂t_sed = 0 - ᶜinv_ρ̂ = (@. lazy( - specific( - FT(1), - Y.c.sgsʲs.:($$j).ρa, - FT(0), - ᶜρʲs.:($$j), - turbconv_model, - ), - )) + ᶜinv_ρ̂ = + @. lazy(specific(FT(1), Y.c.sgsʲs.:($$j).ρa, FT(0), ᶜρʲ, turbconv_model)) # Sedimentation # TODO - lazify ᶜwₗʲs computation. No need to cache it. @@ -462,21 +443,11 @@ function edmfx_sgs_vertical_advection_tendency!( ᶜwʲ = MatrixFields.get_field(p.precomputed, wʲ_name) # Advective form advection of tracers with updraft velocity - va = vertical_advection( - ᶠu³ʲs.:($j), - ᶜqʲ, - edmfx_tracer_upwinding, - ) + va = vertical_advection(ᶠu³ʲs.:($j), ᶜqʲ, edmfx_tracer_upwinding) @. ᶜqʲₜ += va # Flux form sedimentation of tracers - vtt = updraft_sedimentation( - ᶜρʲs.:($j), - ᶜwʲ, - ᶜa, - ᶜqʲ, - ᶠJ, - ) + vtt = updraft_sedimentation(ᶜρʲs.:($j), ᶜwʲ, ᶜa, ᶜqʲ, ᶠJ) @. ᶜqʲₜ += ᶜinv_ρ̂ * vtt @. Yₜ.c.sgsʲs.:($$j).q_tot += ᶜinv_ρ̂ * vtt @. ᶜ∂ρ∂t_sed += vtt @@ -493,15 +464,8 @@ function edmfx_sgs_vertical_advection_tendency!( else error("Unsupported moisture tracer variable") end - vtt = updraft_sedimentation( - ᶜρʲs.:($j), - ᶜwʲ, - ᶜa, - ᶜqʲ .* ᶜmse_li, - ᶠJ, - ) - @. Yₜ.c.sgsʲs.:($$j).mse += - ᶜinv_ρ̂ * vtt + vtt = updraft_sedimentation(ᶜρʲ, ᶜwʲ, ᶜa, ᶜqʲ .* ᶜmse_li, ᶠJ) + @. Yₜ.c.sgsʲs.:($$j).mse += ᶜinv_ρ̂ * vtt end # Contribution of density variation due to sedimentation @@ -542,21 +506,11 @@ function edmfx_sgs_vertical_advection_tendency!( ᶜwʲ = MatrixFields.get_field(p.precomputed, wʲ_name) # Advective form advection of tracers with updraft velocity - va = vertical_transport( - ᶠu³ʲs.:($j), - ᶜχʲ, - edmfx_tracer_upwinding, - ) + va = vertical_transport(ᶠu³ʲs.:($j), ᶜχʲ, edmfx_tracer_upwinding) @. ᶜχʲₜ += va # Flux form sedimentation of tracers - vtt = updraft_sedimentation( - ᶜρʲs.:($j), - ᶜwʲ, - ᶜa, - ᶜχʲ, - ᶠJ, - ) + vtt = updraft_sedimentation(ᶜρʲ, ᶜwʲ, ᶜa, ᶜχʲ, ᶠJ) @. ᶜχʲₜ += ᶜinv_ρ̂ * vtt # Contribution of density variation due to sedimentation diff --git a/src/prognostic_equations/implicit/implicit_tendency.jl b/src/prognostic_equations/implicit/implicit_tendency.jl index 69365901e1..4493316451 100644 --- a/src/prognostic_equations/implicit/implicit_tendency.jl +++ b/src/prognostic_equations/implicit/implicit_tendency.jl @@ -73,28 +73,22 @@ end # the implicit tendency function. Since dt >= dtγ, we can safely use dt for now. function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:none}) - ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J - ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J - return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠu³ * ᶠinterp(ᶜχ)))) + ᶠρ = face_density(ᶜρ) + return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠu³ * ᶠinterp(ᶜχ)))) end function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order}) - ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J - ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J - return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind1(ᶠu³, ᶜχ)))) + ᶠρ = face_density(ᶜρ) + return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠupwind1(ᶠu³, ᶜχ)))) end @static if pkgversion(ClimaCore) ≥ v"0.14.22" function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:vanleer_limiter}) - ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J - ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J - return @. lazy( - -(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠlin_vanleer(ᶠu³, ᶜχ, dt))), - ) + ᶠρ = face_density(ᶜρ) + return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠlin_vanleer(ᶠu³, ᶜχ, dt)))) end end function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order}) - ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J - ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J - return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind3(ᶠu³, ᶜχ)))) + ᶠρ = face_density(ᶜρ) + return @. lazy(-(ᶜadvdivᵥ(ᶠρ * ᶠupwind3(ᶠu³, ᶜχ)))) end vertical_advection(ᶠu³, ᶜχ, ::Val{:none}) = @@ -105,32 +99,26 @@ vertical_advection(ᶠu³, ᶜχ, ::Val{:third_order}) = @. lazy(-(ᶜadvdivᵥ(ᶠupwind3(ᶠu³, ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³))) function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) - (; moisture_model, turbconv_model, rayleigh_sponge, microphysics_model) = - p.atmos + (; moisture_model, turbconv_model, rayleigh_sponge, microphysics_model) = p.atmos (; params, dt) = p n = n_mass_flux_subdomains(turbconv_model) - ᶜJ = Fields.local_geometry_field(axes(Y.c)).J - ᶠJ = Fields.local_geometry_field(axes(Y.f)).J + ᶜρ = Y.c.ρ + ᶠρ = face_density(ᶜρ) (; ᶠgradᵥ_ᶜΦ) = p.core (; ᶠu³, ᶜp, ᶜts) = p.precomputed thermo_params = CAP.thermodynamics_params(params) cp_d = CAP.cp_d(params) - ᶜh_tot = @. lazy( - TD.total_specific_enthalpy( - thermo_params, - ᶜts, - specific(Y.c.ρe_tot, Y.c.ρ), - ), - ) + e_tot = @. lazy(specific(Y.c.ρe_tot, ᶜρ)) + ᶜh_tot = @. lazy(TD.total_specific_enthalpy(thermo_params, ᶜts, e_tot)) - @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠu³) + @. Yₜ.c.ρ -= ᶜdivᵥ(ᶠρ * ᶠu³) # Central vertical advection of active tracers (e_tot and q_tot) - vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜh_tot, dt, Val(:none)) + vtt = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, dt, Val(:none)) @. Yₜ.c.ρe_tot += vtt if !(moisture_model isa DryModel) - ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ)) - vtt = vertical_transport(Y.c.ρ, ᶠu³, ᶜq_tot, dt, Val(:none)) + ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, ᶜρ)) + vtt = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, dt, Val(:none)) @. Yₜ.c.ρq_tot += vtt end @@ -140,63 +128,41 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) # using downward biasing and free outflow bottom boundary condition if moisture_model isa NonEquilMoistModel (; ᶜwₗ, ᶜwᵢ) = p.precomputed - @. Yₜ.c.ρq_liq -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwₗ)) * specific(Y.c.ρq_liq, Y.c.ρ), - ), - ) - @. Yₜ.c.ρq_ice -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwᵢ)) * specific(Y.c.ρq_ice, Y.c.ρ), - ), - ) + q_liq = @. lazy(specific(Y.c.ρq_liq, ᶜρ)) + q_ice = @. lazy(specific(Y.c.ρq_ice, ᶜρ)) + @. Yₜ.c.ρq_liq -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₗ) * q_liq)) + @. Yₜ.c.ρq_ice -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵢ) * q_ice)) end if microphysics_model isa Microphysics1Moment (; ᶜwᵣ, ᶜwₛ) = p.precomputed - @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwᵣ)) * specific(Y.c.ρq_rai, Y.c.ρ), - ), - ) - @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwₛ)) * specific(Y.c.ρq_sno, Y.c.ρ), - ), - ) + q_rai = @. lazy(specific(Y.c.ρq_rai, ᶜρ)) + q_sno = @. lazy(specific(Y.c.ρq_sno, ᶜρ)) + @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵣ) * q_rai)) + @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₛ) * q_sno)) end if microphysics_model isa Microphysics2Moment (; ᶜwₙₗ, ᶜwₙᵣ, ᶜwᵣ, ᶜwₛ) = p.precomputed - @. Yₜ.c.ρn_liq -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwₙₗ)) * specific(Y.c.ρn_liq, Y.c.ρ), - ), - ) - @. Yₜ.c.ρn_rai -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwₙᵣ)) * specific(Y.c.ρn_rai, Y.c.ρ), - ), - ) - @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwᵣ)) * specific(Y.c.ρq_rai, Y.c.ρ), - ), - ) - @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ( - ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias( - Geometry.WVector(-(ᶜwₛ)) * specific(Y.c.ρq_sno, Y.c.ρ), - ), - ) + n_liq = @. lazy(specific(Y.c.ρn_liq, ᶜρ)) + n_rai = @. lazy(specific(Y.c.ρn_rai, ᶜρ)) + q_rai = @. lazy(specific(Y.c.ρq_rai, ᶜρ)) + q_sno = @. lazy(specific(Y.c.ρq_sno, ᶜρ)) + @. Yₜ.c.ρn_liq -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₙₗ) * n_liq)) + @. Yₜ.c.ρn_rai -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₙᵣ) * n_rai)) + @. Yₜ.c.ρq_rai -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵣ) * q_rai)) + @. Yₜ.c.ρq_sno -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwₛ) * q_sno)) end if microphysics_model isa Microphysics2MomentP3 - (; ρ, ρn_ice, ρq_rim, ρb_rim) = Y.c - ᶜwnᵢ = @. lazy(Geometry.WVector(p.precomputed.ᶜwnᵢ)) - ᶜwᵢ = @. lazy(Geometry.WVector(p.precomputed.ᶜwᵢ)) - ᶠρ = @. lazy(ᶠinterp(ρ * ᶜJ) / ᶠJ) + (; ᶜwnᵢ, ᶜwᵢ) = p.precomputed + ᶜρ = Y.c.ρ + ᶠρ = face_density(ᶜρ) # Note: `ρq_ice` is handled above, in `moisture_model isa NonEquilMoistModel` - @. Yₜ.c.ρn_ice -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(- ᶜwnᵢ * specific(ρn_ice, ρ))) - @. Yₜ.c.ρq_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(- ᶜwᵢ * specific(ρq_rim, ρ))) - @. Yₜ.c.ρb_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(- ᶜwᵢ * specific(ρb_rim, ρ))) + n_ice = @. lazy(specific(Y.c.ρn_ice, ᶜρ)) + q_rim = @. lazy(specific(Y.c.ρq_rim, ᶜρ)) + b_rim = @. lazy(specific(Y.c.ρb_rim, ᶜρ)) + @. Yₜ.c.ρn_ice -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwnᵢ) * n_ice)) + @. Yₜ.c.ρq_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵢ) * q_rim)) + @. Yₜ.c.ρb_rim -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-WVec(ᶜwᵢ) * b_rim)) end # TODO - decide if this needs to be explicit or implicit @@ -207,8 +173,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) ᶜθ_v = @. lazy(theta_v(thermo_params, ᶜts)) ᶜθ_vr = @. lazy(theta_vr(thermo_params, ᶜts)) ᶜΠ = @. lazy(dry_exner_function(thermo_params, ᶜts)) - @. Yₜ.f.u₃ -= ᶠgradᵥ_ᶜΦ - ᶠgradᵥ(ᶜΦ_r) + - cp_d * (ᶠinterp(ᶜθ_v - ᶜθ_vr)) * ᶠgradᵥ(ᶜΠ) + @. Yₜ.f.u₃ -= ᶠgradᵥ_ᶜΦ - ᶠgradᵥ(ᶜΦ_r) + cp_d * (ᶠinterp(ᶜθ_v - ᶜθ_vr)) * ᶠgradᵥ(ᶜΠ) if rayleigh_sponge isa RayleighSponge ᶠz = Fields.coordinate_field(Y.f).z @@ -217,8 +182,7 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t) @. Yₜ.f.u₃ -= β_rayleigh_w(rs, ᶠz, zmax) * Y.f.u₃ if turbconv_model isa PrognosticEDMFX for j in 1:n - @. Yₜ.f.sgsʲs.:($$j).u₃ -= - β_rayleigh_w(rs, ᶠz, zmax) * Y.f.sgsʲs.:($$j).u₃ + @. Yₜ.f.sgsʲs.:($$j).u₃ -= β_rayleigh_w(rs, ᶠz, zmax) * Y.f.sgsʲs.:($$j).u₃ end end end diff --git a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl index 4dbd47dd8b..e2f9c37e50 100644 --- a/src/prognostic_equations/implicit/manual_sparse_jacobian.jl +++ b/src/prognostic_equations/implicit/manual_sparse_jacobian.jl @@ -419,6 +419,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) ) ᶜρ = Y.c.ρ + ᶠρ = face_density(ᶜρ) ᶜuₕ = Y.c.uₕ ᶠu₃ = Y.f.u₃ ᶜJ = Fields.local_geometry_field(Y.c).J @@ -453,8 +454,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) @. ᶠp_grad_matrix = DiagonalMatrixRow(-1 / ᶠinterp(ᶜρ)) ⋅ ᶠgradᵥ_matrix() - @. ᶜadvection_matrix = - -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) + @. ᶜadvection_matrix = -(ᶜadvdivᵥ_matrix()) ⋅ DiagonalMatrixRow(ᶠρ) if use_derivative(topography_flag) ∂ᶜρ_err_∂ᶜuₕ = matrix[@name(c.ρ), @name(c.uₕ)] @@ -563,7 +563,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) #TODO: tetsing explicit vs implicit #@. ∂ᶜρe_tot_err_∂ᶜρe_tot += # dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ - # DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅ ᶠright_bias_matrix() ⋅ + # DiagonalMatrixRow(ᶠρ ⋅ ᶠright_bias_matrix() ⋅ # DiagonalMatrixRow( # -(1 + ᶜkappa_m) / ᶜρ * ifelse( # ᶜh_tot == 0, @@ -577,7 +577,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) #TODO: tetsing explicit vs implicit #@. ∂ᶜρe_tot_err_∂ᶜρq_tot = # dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ - # DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅ ᶠright_bias_matrix() ⋅ + # DiagonalMatrixRow(ᶠρ ⋅ ᶠright_bias_matrix() ⋅ # DiagonalMatrixRow( # -(ᶜkappa_m) * ∂e_int_∂q_tot / ᶜρ * ifelse( # ᶜh_tot == 0, @@ -591,7 +591,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) #TODO: testing explicit vs implicit #@. ∂ᶜρq_tot_err_∂ᶜρq_tot = # dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅ - # DiagonalMatrixRow(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) ⋅ ᶠright_bias_matrix() ⋅ + # DiagonalMatrixRow(ᶠρ ⋅ ᶠright_bias_matrix() ⋅ # DiagonalMatrixRow( # -1 / ᶜρ * ifelse( # ᶜq_tot == 0, @@ -1180,10 +1180,9 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) if use_derivative(sgs_mass_flux_flag) # Jacobian contributions of updraft massflux to grid-mean ∂ᶜupdraft_mass_flux_∂ᶜscalar = ᶠbidiagonal_matrix_ct3 + ᶠρ = face_density(ᶜρ) @. ∂ᶜupdraft_mass_flux_∂ᶜscalar = - DiagonalMatrixRow( - (ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) * (ᶠu³ʲs.:(1) - ᶠu³), - ) ⋅ ᶠinterp_matrix() ⋅ + DiagonalMatrixRow(ᶠρ * (ᶠu³ʲs.:(1) - ᶠu³)) ⋅ ᶠinterp_matrix() ⋅ DiagonalMatrixRow(Y.c.sgsʲs.:(1).ρa / ᶜρʲs.:(1)) @. p.scratch.ᶜtridiagonal_matrix_scalar = dtγ * ᶜadvdivᵥ_matrix() ⋅ ∂ᶜupdraft_mass_flux_∂ᶜscalar @@ -1192,8 +1191,7 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t) ## grid-mean ρe_tot ᶜkappa_m = p.scratch.ᶜtemp_scalar @. ᶜkappa_m = - TD.gas_constant_air(thermo_params, ᶜts) / - TD.cv_m(thermo_params, ᶜts) + TD.gas_constant_air(thermo_params, ᶜts) / TD.cv_m(thermo_params, ᶜts) ᶜ∂kappa_m∂q_tot = p.scratch.ᶜtemp_scalar_2 @. ᶜ∂kappa_m∂q_tot = diff --git a/src/prognostic_equations/water_advection.jl b/src/prognostic_equations/water_advection.jl index 231a406e16..66269b85c9 100644 --- a/src/prognostic_equations/water_advection.jl +++ b/src/prognostic_equations/water_advection.jl @@ -47,17 +47,13 @@ Modifies: """ function vertical_advection_of_water_tendency!(Yₜ, Y, p, t) - ᶜJ = Fields.local_geometry_field(Y.c).J - ᶠJ = Fields.local_geometry_field(Y.f).J + ᶠρ = face_density(Y.c.ρ) (; ᶜwₜqₜ, ᶜwₕhₜ) = p.precomputed if !(p.atmos.moisture_model isa DryModel) - @. Yₜ.c.ρ -= - ᶜprecipdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(-(ᶜwₜqₜ))) - @. Yₜ.c.ρe_tot -= - ᶜprecipdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(-(ᶜwₕhₜ))) - @. Yₜ.c.ρq_tot -= - ᶜprecipdivᵥ(ᶠinterp(Y.c.ρ * ᶜJ) / ᶠJ * ᶠright_bias(-(ᶜwₜqₜ))) + @. Yₜ.c.ρ -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-(ᶜwₜqₜ))) + @. Yₜ.c.ρe_tot -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-(ᶜwₕhₜ))) + @. Yₜ.c.ρq_tot -= ᶜprecipdivᵥ(ᶠρ * ᶠright_bias(-(ᶜwₜqₜ))) end return nothing diff --git a/src/utils/abbreviations.jl b/src/utils/abbreviations.jl index bedfec0a32..604459a6b3 100644 --- a/src/utils/abbreviations.jl +++ b/src/utils/abbreviations.jl @@ -13,6 +13,8 @@ const CT12 = Geometry.Contravariant12Vector const CT3 = Geometry.Contravariant3Vector const CT123 = Geometry.Contravariant123Vector const UVW = Geometry.UVWVector +const UV = Geometry.UVVector +const WVec = Geometry.WVector const divₕ = Operators.Divergence() const wdivₕ = Operators.WeakDivergence() diff --git a/src/utils/utilities.jl b/src/utils/utilities.jl index fe0b74a10a..38bba65db4 100644 --- a/src/utils/utilities.jl +++ b/src/utils/utilities.jl @@ -192,6 +192,22 @@ function strain_rate_norm(S, axis = Geometry.UVWAxis()) return S_norm end +""" + face_density(ᶜρ) + +Compute the density at cell faces from the density at cell centers. + +This reconstruction satisfies a Jacobian-weighted average, see the definition in +Eq. (104) in Section 5.6 of the dycore paper (Yatunin et al., 2025). +""" +function face_density(ᶜρ) + ᶠspace = Spaces.face_space(axes(ᶜρ)) + ᶠJ = Fields.local_geometry_field(ᶠspace).J + ᶜJ = Fields.local_geometry_field(ᶜρ).J + + @. lazy(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ) +end + """ g³³_field(space)