From e05421837f3e68b543f579f62d97eb35b87b959c Mon Sep 17 00:00:00 2001 From: Zhaoyi Shen <11598433+szy21@users.noreply.github.com> Date: Thu, 12 Oct 2023 21:35:21 -0700 Subject: [PATCH] wip --- .../advective_edmfx_bomex_box.yml | 2 +- src/prognostic_equations/advection.jl | 8 +- src/prognostic_equations/edmfx_closures.jl | 4 +- src/prognostic_equations/hyperdiffusion.jl | 77 +++++++++++++++++-- 4 files changed, 76 insertions(+), 15 deletions(-) diff --git a/config/model_configs/advective_edmfx_bomex_box.yml b/config/model_configs/advective_edmfx_bomex_box.yml index 116c7be8a7c..f53c8b87d5c 100644 --- a/config/model_configs/advective_edmfx_bomex_box.yml +++ b/config/model_configs/advective_edmfx_bomex_box.yml @@ -24,7 +24,7 @@ y_elem: 2 z_elem: 60 z_stretch: false perturb_initstate: false -dt: "2secs" +dt: "5secs" t_end: "6hours" dt_save_to_disk: "1mins" toml: [toml/edmfx_box.toml] diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl index 8f15c1ae7f4..b08c268460f 100644 --- a/src/prognostic_equations/advection.jl +++ b/src/prognostic_equations/advection.jl @@ -44,12 +44,12 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t) end end - # FIXME: add horizontal advection if p.atmos.turbconv_model isa AdvectiveEDMFX for j in 1:n - if :h_tot in propertynames(Y.c) - @. Yₜ.c.sgsʲs.:($$j).h_tot -= 0 - end + @. Yₜ.c.sgsʲs.:($$j).h_tot -= + adjoint(CA.CT12(Y.c.uₕ)) * CA.gradₕ(Y.c.sgsʲs.:($$j).h_tot) + @. Yₜ.c.sgsʲs.:($$j).q_tot -= + adjoint(CA.CT12(Y.c.uₕ)) * CA.gradₕ(Y.c.sgsʲs.:($$j).q_tot) end end diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl index 514a79746cb..ce53a8614fb 100644 --- a/src/prognostic_equations/edmfx_closures.jl +++ b/src/prognostic_equations/edmfx_closures.jl @@ -69,8 +69,8 @@ function ᶠupdraft_nh_pressure( #plume_scale_height = max(ᶠa * 10e3, H_up_min) # We also used to have advection term here: α_a * w_up * div_w_up - return 2 * α_d * (ᶠu3ʲ - ᶠu3⁰) * CC.Geometry._norm(ᶠu3ʲ - ᶠu3⁰, ᶠlg) / - plume_scale_height #+ α_b * ᶠbuoyʲ + return 1.2 * α_d * (ᶠu3ʲ - ᶠu3⁰) * CC.Geometry._norm(ᶠu3ʲ - ᶠu3⁰, ᶠlg) / + plume_scale_height + α_b * ᶠbuoyʲ end end diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl index a323383d09d..e2b40646b9e 100644 --- a/src/prognostic_equations/hyperdiffusion.jl +++ b/src/prognostic_equations/hyperdiffusion.jl @@ -21,7 +21,10 @@ function hyperdiffusion_cache(Y, atmos, do_dss) # Sub-grid scale quantities ᶜ∇²uʲs = - atmos.turbconv_model isa EDMFX ? similar(Y.c, NTuple{n, C123{FT}}) : (;) + ( + atmos.turbconv_model isa EDMFX || + atmos.turbconv_model isa AdvectiveEDMFX + ) ? similar(Y.c, NTuple{n, C123{FT}}) : (;) sgs_quantities = atmos.turbconv_model isa EDMFX ? (; @@ -33,6 +36,14 @@ function hyperdiffusion_cache(Y, atmos, do_dss) specific_sgsʲs.(Y.c, atmos.turbconv_model) ), ) : + atmos.turbconv_model isa AdvectiveEDMFX ? + (; + ᶜ∇²tke⁰ = similar(Y.c, FT), + ᶜ∇²uₕʲs = similar(Y.c, NTuple{n, C12{FT}}), + ᶜ∇²uᵥʲs = similar(Y.c, NTuple{n, C3{FT}}), + ᶜ∇²h_totʲs = similar(Y.c, NTuple{n, FT}), + ᶜ∇²q_totʲs = similar(Y.c, NTuple{n, FT}), + ) : atmos.turbconv_model isa DiagnosticEDMFX ? (; ᶜ∇²tke⁰ = similar(Y.c, FT)) : (;) quantities = (; gs_quantities..., sgs_quantities...) @@ -71,6 +82,9 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) ᶜ∇²specific_energyʲs, ) = p end + if turbconv_model isa AdvectiveEDMFX + (; ᶜρa⁰, ᶜρʲs, ᶜ∇²tke⁰, ᶜtke⁰, ᶜ∇²uₕʲs, ᶜ∇²uᵥʲs, ᶜ∇²uʲs, ᶜ∇²h_totʲs) = p + end if turbconv_model isa DiagnosticEDMFX (; ᶜtke⁰, ᶜ∇²tke⁰) = p end @@ -89,7 +103,9 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) if turbconv_model isa EDMFX && diffuse_tke @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜspecific⁰.tke)) end - if turbconv_model isa DiagnosticEDMFX && diffuse_tke + if ( + turbconv_model isa DiagnosticEDMFX || turbconv_model isa AdvectiveEDMFX + ) && diffuse_tke @. ᶜ∇²tke⁰ = wdivₕ(gradₕ(ᶜtke⁰)) end @@ -110,6 +126,15 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) end end + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + @. ᶜ∇²uʲs.:($$j) = + C123(wgradₕ(divₕ(p.ᶜuʲs.:($$j)))) - + C123(wcurlₕ(C123(curlₕ(p.ᶜuʲs.:($$j))))) + @. ᶜ∇²h_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).h_tot)) + end + end + if do_dss NVTX.@range "dss_hyperdiffusion_tendency" color = colorant"green" begin for dss_op! in ( @@ -122,10 +147,11 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) # operations do not accept Covariant123Vector types dss_op!(ᶜ∇²u, buffer.ᶜ∇²u) dss_op!(ᶜ∇²specific_energy, buffer.ᶜ∇²specific_energy) - if turbconv_model isa EDMFX && diffuse_tke + if diffuse_tke dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰) end - if turbconv_model isa EDMFX + if (turbconv_model isa EDMFX) || + (turbconv_model isa AdvectiveEDMFX) # Need to split the DSS computation here, because our DSS # operations do not accept Covariant123Vector types for j in 1:n @@ -138,11 +164,16 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) @. ᶜ∇²uʲs.:($$j) = C123(ᶜ∇²uₕʲs.:($$j)) + C123(ᶜ∇²uᵥʲs.:($$j)) end + end + if turbconv_model isa EDMFX dss_op!(ᶜ∇²specific_energyʲs, buffer.ᶜ∇²specific_energyʲs) end - if turbconv_model isa DiagnosticEDMFX && diffuse_tke - dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰) + if turbconv_model isa AdvectiveEDMFX + dss_op!(ᶜ∇²h_totʲs, buffer.ᶜ∇²h_totʲs) end + # if turbconv_model isa DiagnosticEDMFX && diffuse_tke + # dss_op!(ᶜ∇²tke⁰, buffer.ᶜ∇²tke⁰) + # end end end end @@ -158,7 +189,8 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) @. ᶜρ_energyₜ -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²specific_energy)) # Sub-grid scale hyperdiffusion continued - if turbconv_model isa EDMFX && diffuse_tke + if (turbconv_model isa EDMFX || turbconv_model isa AdvectiveEDMFX) && + diffuse_tke @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(ᶜρa⁰ * gradₕ(ᶜ∇²tke⁰)) end if turbconv_model isa EDMFX @@ -177,7 +209,17 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t) wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²specific_energyʲs.:($$j))) end end - # FIXME: add hyperdiffusion for advective edmf + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + if point_type <: Geometry.Abstract3DPoint + # only need curl-curl part + @. ᶜ∇²uᵥʲs.:($$j) = C3(wcurlₕ(C123(curlₕ(ᶜ∇²uʲs.:($$j))))) + @. Yₜ.f.sgsʲs.:($$j).u₃ += + κ₄ * ᶠwinterp(ᶜJ * Y.c.ρ, ᶜ∇²uᵥʲs.:($$j)) + end + @. Yₜ.c.sgsʲs.:($$j).h_tot -= κ₄ * wdivₕ(gradₕ(ᶜ∇²h_totʲs.:($$j))) + end + end if turbconv_model isa DiagnosticEDMFX && diffuse_tke @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²tke⁰)) @@ -195,6 +237,9 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) if turbconv_model isa EDMFX (; ᶜspecificʲs, ᶜ∇²specific_tracersʲs, ᶜp) = p end + if turbconv_model isa EDMFX + (; ᶜq_totʲs) = p + end if do_dss buffer = p.hyperdiffusion_ghost_buffer end @@ -211,6 +256,12 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) end end + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + @. ᶜ∇²q_totʲs.:($$j) = wdivₕ(gradₕ(Y.c.sgsʲs.:($$j).q_tot)) + end + end + if do_dss NVTX.@range "dss_hyperdiffusion_tendency" color = colorant"green" begin for dss_op! in ( @@ -222,6 +273,9 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) if turbconv_model isa EDMFX dss_op!(ᶜ∇²specific_tracersʲs, buffer.ᶜ∇²specific_tracersʲs) end + if turbconv_model isa EDMFX + dss_op!(ᶜ∇²q_totʲs, buffer.ᶜ∇²q_totʲs) + end end end end @@ -247,5 +301,12 @@ NVTX.@annotate function tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t) end end end + if turbconv_model isa AdvectiveEDMFX + for j in 1:n + @. Yₜ.c.sgsʲs.:($$j).ρa -= + κ₄ * wdivₕ(Y.c.sgsʲs.:($$j).ρa * gradₕ(ᶜ∇²q_totʲs.:($$j))) + @. Yₜ.c.sgsʲs.:($$j).q_tot -= κ₄ * wdivₕ(gradₕ(ᶜ∇²q_totʲs.:($$j))) + end + end return nothing end