wip

CliMA · Oct 10, 2023 · f8853c4 · f8853c4
1 parent 8384a63
commit f8853c4
Show file tree

Hide file tree

Showing 11 changed files with 74 additions and 47 deletions.
diff --git a/src/cache/advective_edmf_precomputed_quantities.jl b/src/cache/advective_edmf_precomputed_quantities.jl
@@ -25,15 +25,7 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     (; ᶜspecific, ᶜp, ᶜΦ, ᶜh_tot, ᶜρ_ref) = p
     (; ᶜtke⁰, ᶜρa⁰, ᶠu₃⁰, ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶜts⁰, ᶜρ⁰, ᶜh_tot⁰, ᶜq_tot⁰) = p
     (; ᶜmixing_length, ᶜlinear_buoygrad, ᶜstrain_rate_norm, ᶜK_u, ᶜK_h) = p
-    (;
-        ᶜuʲs,
-        ᶠu³ʲs,
-        ᶜKʲs,
-        ᶜtsʲs,
-        ᶜρʲs,
-        ᶜentrʲs,
-        ᶜdetrʲs,
-    ) = p
+    (; ᶜuʲs, ᶠu³ʲs, ᶜKʲs, ᶜtsʲs, ᶜρʲs, ᶜentrʲs, ᶜdetrʲs) = p
     (; ustar, obukhov_length, buoyancy_flux) = p.sfc_conditions
 
     @. ᶜρa⁰ = ρa⁰(Y.c)
@@ -56,12 +48,7 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
     set_sgs_ᶠu₃!(u₃⁰, ᶠu₃⁰, Y, turbconv_model)
     set_velocity_quantities!(ᶜu⁰, ᶠu³⁰, ᶜK⁰, ᶠu₃⁰, Y.c.uₕ, ᶠuₕ³)
     @. ᶜK⁰ += ᶜtke⁰
-    @. ᶜts⁰ = TD.PhaseEquil_phq(
-        thermo_params,
-        ᶜp,
-        ᶜh_tot⁰ - ᶜK⁰ - ᶜΦ,
-        ᶜq_tot⁰,
-    )
+    @. ᶜts⁰ = TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_tot⁰ - ᶜK⁰ - ᶜΦ, ᶜq_tot⁰)
     @. ᶜρ⁰ = TD.air_density(thermo_params, ᶜts⁰)
 
     for j in 1:n
@@ -75,13 +62,8 @@ function set_advective_edmf_precomputed_quantities!(Y, p, ᶠuₕ³, t)
         ᶜq_totʲ = Y.c.sgsʲs.:($j).q_tot
 
         set_velocity_quantities!(ᶜuʲ, ᶠu³ʲ, ᶜKʲ, ᶠu₃ʲ, Y.c.uₕ, ᶠuₕ³)
-        # FIXME: get ts
-        @. ᶜtsʲ = TD.PhaseEquil_phq(
-            thermo_params,
-            ᶜp,
-            ᶜh_totʲ - ᶜKʲ - ᶜΦ,
-            ᶜq_totʲ,
-        )
+        @. ᶜtsʲ =
+            TD.PhaseEquil_phq(thermo_params, ᶜp, ᶜh_totʲ - ᶜKʲ - ᶜΦ, ᶜq_totʲ)
         @. ᶜρʲ = TD.air_density(thermo_params, ᶜtsʲ)
 
         # When ᶜe_intʲ = ᶜe_int and ᶜq_totʲ = ᶜq_tot, we still observe that

diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
@@ -147,7 +147,12 @@ function precomputed_quantities(Y, atmos)
             ᶜK_h = similar(Y.c, FT),
             ρatke_flux = similar(Fields.level(Y.f, half), C3{FT}),
         ) : (;)
-    return (; gs_quantities..., sgs_quantities..., advective_sgs_quantities..., diagnostic_sgs_quantities...)
+    return (;
+        gs_quantities...,
+        sgs_quantities...,
+        advective_sgs_quantities...,
+        diagnostic_sgs_quantities...,
+    )
 end
 
 # Interpolates the third contravariant component of Y.c.uₕ to cell faces.

diff --git a/src/initial_conditions/atmos_state.jl b/src/initial_conditions/atmos_state.jl
@@ -144,9 +144,9 @@ function turbconv_center_variables(ls, turbconv_model::AdvectiveEDMFX, gs_vars)
     sgs⁰ = (; ρatke = ls.ρ * (1 - a_draft) * ls.turbconv_state.tke)
     ρa = ls.ρ * a_draft / n
     h_tot =
-            TD.specific_enthalpy(ls.thermo_params, ls.thermo_state) +
-            norm_sqr(ls.velocity) / 2 +
-            CAP.grav(ls.params) * ls.geometry.coordinates.z
+        TD.specific_enthalpy(ls.thermo_params, ls.thermo_state) +
+        norm_sqr(ls.velocity) / 2 +
+        CAP.grav(ls.params) * ls.geometry.coordinates.z
     q_tot = TD.total_specific_humidity(ls.thermo_params, ls.thermo_state)
     sgsʲs = ntuple(_ -> (; ρa = ρa, h_tot = h_tot, q_tot = q_tot), Val(n))
     return (; sgs⁰, sgsʲs)

diff --git a/src/prognostic_equations/advection.jl b/src/prognostic_equations/advection.jl
@@ -12,12 +12,14 @@ NVTX.@annotate function horizontal_advection_tendency!(Yₜ, Y, p, t)
     if p.atmos.turbconv_model isa AbstractEDMF
         (; ᶜu⁰) = p
     end
-    if p.atmos.turbconv_model isa EDMFX || p.atmos.turbconv_model isa AdvectiveEDMFX
+    if p.atmos.turbconv_model isa EDMFX ||
+       p.atmos.turbconv_model isa AdvectiveEDMFX
         (; ᶜuʲs) = p
     end
 
     @. Yₜ.c.ρ -= wdivₕ(Y.c.ρ * ᶜu)
-    if p.atmos.turbconv_model isa EDMFX || p.atmos.turbconv_model isa AdvectiveEDMFX
+    if p.atmos.turbconv_model isa EDMFX ||
+       p.atmos.turbconv_model isa AdvectiveEDMFX
         for j in 1:n
             @. Yₜ.c.sgsʲs.:($$j).ρa -= wdivₕ(Y.c.sgsʲs.:($$j).ρa * ᶜuʲs.:($$j))
         end
@@ -108,7 +110,9 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     (; ᶠu³⁰) = advect_tke ? p : all_nothing
     ᶜρa⁰ = advect_tke ? (n > 0 ? p.ᶜρa⁰ : Y.c.ρ) : nothing
     ᶜρ⁰ = advect_tke ? (n > 0 ? p.ᶜρ⁰ : Y.c.ρ) : nothing
-    ᶜtke⁰ = advect_tke ? (turbconv_model isa EDMFX ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) : nothing
+    ᶜtke⁰ =
+        advect_tke ? (turbconv_model isa EDMFX ? p.ᶜspecific⁰.tke : p.ᶜtke⁰) :
+        nothing
     ᶜa_scalar = p.ᶜtemp_scalar
     ᶜω³ = p.ᶜtemp_CT3
     ᶠω¹² = p.ᶠtemp_CT12
@@ -173,8 +177,10 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
 
                 if :ρae_tot in propertynames(Yₜ.c.sgsʲs.:($j))
                     @. ᶜa_scalar[colidx] =
-                        ᶜh_totʲs.:($$j)[colidx] *
-                        draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx])
+                        ᶜh_totʲs.:($$j)[colidx] * draft_area(
+                            Y.c.sgsʲs.:($$j).ρa[colidx],
+                            ᶜρʲs.:($$j)[colidx],
+                        )
                     vertical_transport!(
                         Yₜ.c.sgsʲs.:($j).ρae_tot[colidx],
                         ᶜJ[colidx],
@@ -190,8 +196,10 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
                     matching_subfields(Yₜ.c.sgsʲs.:($j), ᶜspecificʲs.:($j))
                     χ_name == :e_tot && continue
                     @. ᶜa_scalar[colidx] =
-                        ᶜχʲ[colidx] *
-                        draft_area(Y.c.sgsʲs.:($$j).ρa[colidx], ᶜρʲs.:($$j)[colidx])
+                        ᶜχʲ[colidx] * draft_area(
+                            Y.c.sgsʲs.:($$j).ρa[colidx],
+                            ᶜρʲs.:($$j)[colidx],
+                        )
                     vertical_transport!(
                         ᶜρaχʲₜ[colidx],
                         ᶜJ[colidx],
@@ -220,11 +228,15 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
                     edmfx_upwinding,
                 )
 
-                if :h_tot in propertynames(Yₜ.c.sgsʲs.:($j))
-                    @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] += 0
-                end
-
-                @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] += 0
+                # FIXME: boundary conditions
+                @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] -= ᶜinterp(
+                    adjoint(ᶠu³ʲs.:($$j)[colidx]) *
+                    ᶠgradᵥ(Y.c.sgsʲs.:($$j).h_tot[colidx]),
+                )
+                @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] -= ᶜinterp(
+                    adjoint(ᶠu³ʲs.:($$j)[colidx]) *
+                    ᶠgradᵥ(Y.c.sgsʲs.:($$j).q_tot[colidx]),
+                )
             end
         end
 

diff --git a/src/prognostic_equations/edmfx_entr_detr.jl b/src/prognostic_equations/edmfx_entr_detr.jl
@@ -315,7 +315,14 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX)
     return nothing
 end
 
-function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::AdvectiveEDMFX)
+function edmfx_entr_detr_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    turbconv_model::AdvectiveEDMFX,
+)
 
     n = n_mass_flux_subdomains(turbconv_model)
     (; ᶜentrʲs, ᶜdetrʲs) = p
@@ -330,7 +337,7 @@ function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::Advect
         @. Yₜ.c.sgsʲs.:($$j).h_tot[colidx] +=
             ᶜentrʲs.:($$j)[colidx] *
             (ᶜh_tot⁰[colidx] - Y.c.sgsʲs.:($$j).h_tot[colidx])
-        
+
         @. Yₜ.c.sgsʲs.:($$j).q_tot[colidx] +=
             ᶜentrʲs.:($$j)[colidx] *
             (ᶜq_tot⁰[colidx] - Y.c.sgsʲs.:($$j).q_tot[colidx])

diff --git a/src/prognostic_equations/edmfx_sgs_flux.jl b/src/prognostic_equations/edmfx_sgs_flux.jl
@@ -90,8 +90,14 @@ function edmfx_sgs_mass_flux_tendency!(
 end
 
 # FIXME: add sgs mass flux
-edmfx_sgs_mass_flux_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::AdvectiveEDMFX) =
-    nothing
+edmfx_sgs_mass_flux_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    turbconv_model::AdvectiveEDMFX,
+) = nothing
 
 function edmfx_sgs_mass_flux_tendency!(
     Yₜ,
@@ -225,7 +231,14 @@ function edmfx_sgs_diffusive_flux_tendency!(
 end
 
 # FIXME: add SGS diffusive flux
-edmfx_sgs_diffusive_flux_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::AdvectiveEDMFX) = nothing
+edmfx_sgs_diffusive_flux_tendency!(
+    Yₜ,
+    Y,
+    p,
+    t,
+    colidx,
+    turbconv_model::AdvectiveEDMFX,
+) = nothing
 
 function edmfx_sgs_diffusive_flux_tendency!(
     Yₜ,

diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl
@@ -178,7 +178,7 @@ NVTX.@annotate function hyperdiffusion_tendency!(Yₜ, Y, p, t)
         end
     end
     # FIXME: add hyperdiffusion for advective edmf
-    
+
     if turbconv_model isa DiagnosticEDMFX && diffuse_tke
         @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²tke⁰))
     end

diff --git a/src/prognostic_equations/zero_velocity.jl b/src/prognostic_equations/zero_velocity.jl
@@ -10,7 +10,8 @@ function zero_velocity_tendency!(Yₜ, Y, p, t, colidx)
 
         @. Yₜ.c.uₕ[colidx] = C12(FT(0), FT(0))
         @. Yₜ.f.u₃[colidx] = Geometry.Covariant3Vector(FT(0))
-        if p.atmos.turbconv_model isa EDMFX || p.atmos.turbconv_model isa AdvectiveEDMFX 
+        if p.atmos.turbconv_model isa EDMFX ||
+           p.atmos.turbconv_model isa AdvectiveEDMFX
             for j in 1:n
                 @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] =
                     Geometry.Covariant3Vector(FT(0))

diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl
@@ -341,7 +341,8 @@ function get_turbconv_model(
     turbconv_params,
 )
     turbconv = parsed_args["turbconv"]
-    @assert turbconv in (nothing, "edmf", "edmfx", "advective_edmfx", "diagnostic_edmfx")
+    @assert turbconv in
+            (nothing, "edmf", "edmfx", "advective_edmfx", "diagnostic_edmfx")
 
     return if turbconv == "edmf"
         TC.EDMFModel(

diff --git a/src/solver/types.jl b/src/solver/types.jl
@@ -160,7 +160,8 @@ EDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = EDMFX{N, TKE, FT}(a_half)
 struct AdvectiveEDMFX{N, TKE, FT} <: AbstractEDMF
     a_half::FT # WARNING: this should never be used outside of divide_by_ρa
 end
-AdvectiveEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} = AdvectiveEDMFX{N, TKE, FT}(a_half)
+AdvectiveEDMFX{N, TKE}(a_half::FT) where {N, TKE, FT} =
+    AdvectiveEDMFX{N, TKE, FT}(a_half)
 
 struct DiagnosticEDMFX{N, TKE, FT} <: AbstractEDMF
     a_int::FT # area fraction of the first interior cell above the surface

diff --git a/toml/edmfx_box.toml b/toml/edmfx_box.toml
@@ -13,3 +13,8 @@ alias = "min_area"
 value = 1.0e-3
 type = "float"
 description = "Minimum area fraction per updraft. Parameter not described in the literature."
+
+[entr_coeff]
+alias = "entr_coeff"
+value = 1.0e-1
+type = "float"
-Original file line number
+Diff line change
@@ Expand Up @@
             end
         end
         # FIXME: add hyperdiffusion for advective edmf
         if turbconv_model isa DiagnosticEDMFX && diffuse_tke
             @. Yₜ.c.sgs⁰.ρatke -= κ₄ * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²tke⁰))
         end
@@ Expand Down @@