make it work

Author: haakon-e

config/model_configs/kinematic_driver.yml

@@ -2,6 +2,7 @@
 prescribed_flow: true
 initial_condition: "ShipwayHill2012"
 surface_setup: "ShipwayHill2012"
+energy_q_tot_upwinding: first_order_kid
 tracer_upwinding: first_order
 # moist: "nonequil"
 # precip_model: "1M"
@@ -12,13 +13,14 @@ call_cloud_diagnostics_per_stage: true
 config: "column"
 hyperdiff: "false"
 ## Simulation
-z_max: 4e3
-z_elem: 200
+z_max: 2e3
+z_elem: 128
 z_stretch: false
+use_auto_jacobian: true # Needed!! (+only 1 Newton iteration, which is the default)
 # ode_algo: "ARS111" # try: Explicit Euler
 dt: "1secs"
-# t_end: "1hours"
-t_end: "20mins"
+t_end: "1hours"
+# t_end: "20mins"
 dt_save_state_to_disk: "1hours"
 toml: [toml/kinematic_driver.toml]
 check_nan_every: 1
@@ -27,7 +29,7 @@ output_default_diagnostics: false
 diagnostics:
   - short_name: [
       # (z,t) fields
-      ta, thetaa, ha, rhoa, wa, hur, hus, cl, clw, cli,
+      ta, thetaa, ha, rhoa, wa, hur, hus, cl, clw, cli, ua, va, ke,
       # (t,) fields
       lwp, pr,
       # 1M

post_processing/ci_plots.jl

@@ -1699,7 +1699,7 @@ function make_plots(::Val{:kinematic_driver}, output_paths::Vector{<:AbstractStr
         return var
     end
     simdirs = SimDir.(output_paths)
-    short_names = ["hus", "clw", "husra", "ta", "thetaa", "wa", "cli", "hussn"]
+    short_names = ["hus", "clw", "husra", "rhoa", "ta", "thetaa", "wa", "cli", "hussn", "ua", "va", "ke"]
     short_names = short_names ∩ collect(keys(simdirs[1].vars))
     vars = map_comparison(simdirs, short_names) do simdir, short_name
         var = slice(get(simdir; short_name), x = 0, y = 0)

src/cache/precomputed_quantities.jl

@@ -480,8 +480,10 @@ NVTX.@annotate function set_implicit_precomputed_quantities_part1!(Y, p, t)
 
     # TODO: We might want to move this to dss! (and rename dss! to something
     # like enforce_constraints!).
-    set_velocity_at_surface!(Y, ᶠuₕ³, turbconv_model)
-    set_velocity_at_top!(Y, turbconv_model)
+    if !(p.atmos.prescribed_flow isa PrescribedFlow)
+        set_velocity_at_surface!(Y, ᶠuₕ³, turbconv_model)
+        set_velocity_at_top!(Y, turbconv_model)
+    end
 
     set_velocity_quantities!(ᶜu, ᶠu³, ᶜK, Y.f.u₃, Y.c.uₕ, ᶠuₕ³)
     ᶜJ = Fields.local_geometry_field(Y.c).J

src/prognostic_equations/advection.jl

@@ -286,6 +286,10 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
         vtt = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, FT(dt), energy_q_tot_upwinding)
         vtt_central = vertical_transport(ᶜρ, ᶠu³, ᶜq_tot, FT(dt), Val(:none))
+        if energy_q_tot_upwinding == Val(:first_order_kid)
+            vtt = vertical_transport_kid(ᶜρ, ᶠu³, ᶜq_tot, FT(dt), FT(t))
+            vtt_central = NullBroadcasted()  # turn off implicit advection
+        end
         @. Yₜ.c.ρq_tot += vtt - vtt_central
     end
 

src/prognostic_equations/dss.jl

@@ -87,26 +87,25 @@ end
 
 prescribe_flow!(Y, p, t, ::Nothing) = nothing
 function prescribe_flow!(Y, p, t, flow::PrescribedFlow)
-    FT = eltype(p.params)
-    (; prescribed_u₃) = flow
     ᶠlg = Fields.local_geometry_field(Y.f)
-    @. Y.f.u₃ = C3(Geometry.WVector(prescribed_u₃(FT, t)), ᶠlg)
+    z = Fields.coordinate_field(Y.f).z
+    @. Y.f.u₃ = C3(Geometry.WVector(flow(z, t)), ᶠlg)
 
-    return nothing  # comment out to try fixing energy
+    # return nothing  # comment out to try fixing energy
 
     ### Fix energy to initial temperature
     ᶜlg = Fields.local_geometry_field(Y.c)
     local_state = InitialConditions.ShipwayHill2012()(p.params)
-    get_ρ_init(ls) = TD.air_density(ls.thermo_params, ls.thermo_state)
+    get_ρ_init_dry(ls) = ls.thermo_state.ρ * (1 - ls.thermo_state.q_tot)
     get_T_init(ls) = TD.air_temperature(ls.thermo_params, ls.thermo_state)
-    ᶜρ_init = @. lazy(get_ρ_init(local_state(ᶜlg)))
+    ᶜρ_init_dry = @. lazy(get_ρ_init_dry(local_state(ᶜlg)))
     ᶜT_init = @. lazy(get_T_init(local_state(ᶜlg)))
 
     thermo_params = CAP.thermodynamics_params(p.params)
     grav = CAP.grav(p.params)
     z = Fields.coordinate_field(Y.c).z
 
-    @. Y.c.ρ = ᶜρ_init + Y.c.ρq_tot
+    @. Y.c.ρ = ᶜρ_init_dry + Y.c.ρq_tot
     ᶜts = @. lazy(TD.PhaseEquil_ρTq(thermo_params, Y.c.ρ, ᶜT_init, Y.c.ρq_tot / Y.c.ρ))
     ᶜe_kin = compute_kinetic(Y.c.uₕ, Y.f.u₃)
     ᶜe_pot = @. lazy(grav * z)

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -82,6 +82,22 @@ function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order})
     ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
     return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind1(ᶠu³, ᶜχ))))
 end
+vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:first_order_kid}) =
+    vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, Val(:first_order))  # transport e.g. ρe_tot as usual
+function vertical_transport_kid(ᶜρ, ᶠu³, ᶜχ, dt, t)
+    FT = eltype(ᶜρ)
+    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
+    ᶠJ = Fields.local_geometry_field(axes(ᶠu³)).J
+    ρ_sfc = FT(1.1650101)  # TODO: Get from initial conditions or state
+    q_tot_sfc = FT(0.014778325)  # TODO: Get from initial conditions or state
+    u₃_sfc = ShipwayHill2012VelocityProfile{FT}()(FT(0), t)
+    
+    bottom_bc = Geometry.WVector(ρ_sfc * q_tot_sfc * u₃_sfc)
+    ᶜadvdivᵥ_kid = Operators.DivergenceF2C(
+        bottom = Operators.SetValue(bottom_bc), top = Operators.Extrapolate(),
+    )
+    return @. lazy(-(ᶜadvdivᵥ_kid(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠ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
@@ -131,6 +147,9 @@ function implicit_vertical_advection_tendency!(Yₜ, Y, p, t)
     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))
+        if p.atmos.numerics.energy_q_tot_upwinding == Val(:first_order_kid)
+            vtt = NullBroadcasted()  # Turn off implicit energy q_tot advection
+        end
         @. Yₜ.c.ρq_tot += vtt
     end
 

src/solver/type_getters.jl

@@ -115,13 +115,7 @@ function get_atmos(config::AtmosConfig, params)
     )
 
     if parsed_args["prescribed_flow"]
-        function prescribed_u₃(FT::Type{<:Real}, _t)
-            t = FT(_t)
-            w1 = FT(2)
-            t1 = FT(600)  # 10 minutes
-            return t < t1 ? w1 * sin(π * t / t1) : FT(0)
-        end
-        prescribed_flow = PrescribedFlow(; prescribed_u₃)
+        prescribed_flow = ShipwayHill2012VelocityProfile{FT}()
     else
         prescribed_flow = nothing
     end

src/solver/types.jl

@@ -6,11 +6,6 @@ import Dates
 import ClimaUtilities.ClimaArtifacts: @clima_artifact
 import LazyArtifacts
 
-abstract type AbstractPrescribedFlow end
-@kwdef struct PrescribedFlow
-    prescribed_u₃::Function
-end
-
 abstract type AbstractMoistureModel end
 abstract type AbstractMoistModel <: AbstractMoistureModel end
 struct DryModel <: AbstractMoistureModel end
@@ -504,6 +499,15 @@ struct RadiationTRMM_LBA{R}
     end
 end
 
+abstract type PrescribedFlow end
+
+struct ShipwayHill2012VelocityProfile{FT} <: PrescribedFlow end
+function (::ShipwayHill2012VelocityProfile{FT})(z, t) where {FT}
+    w1 = FT(1.5)
+    t1 = FT(600)
+    return t < t1 ? w1 * sinpi(FT(t) / t1) : FT(0)
+end
+
 struct TestDycoreConsistency end
 
 abstract type AbstractTimesteppingMode end
@@ -708,7 +712,7 @@ const ATMOS_MODEL_GROUPS = (
     (AtmosWater, :water),
     (AtmosRadiation, :radiation),
     (AtmosTurbconv, :turbconv),
-    (PrescribedFlow, :prescribed_flow),
+    (ShipwayHill2012VelocityProfile, :prescribed_flow),
     (AtmosGravityWave, :gravity_wave),
     (AtmosSponge, :sponge),
     (AtmosSurface, :surface),

src/surface_conditions/surface_setups.jl

@@ -254,8 +254,10 @@ function (::ShipwayHill2012)(params)
         # Adjust the coefficients from 20 m to the actual value of z.
         z0 = FT(1.5e-4)
         adjustment = (log(20 / z0) / log(interior_z / z0))^2
-        Cd = FT(0.001229) * adjustment
-        Ch = FT(0.001094) * adjustment
+        # Cd = FT(0.001229) * adjustment
+        # Ch = FT(0.001094) * adjustment
+        Cd = FT(0.0)
+        Ch = FT(0.0)
         parameterization = ExchangeCoefficients(; Cd, Ch)
         # <<<
         return SurfaceState(; parameterization, T, p, q_vap)