wip, [ci skip]

CliMA · Jan 29, 2025 · ffff256 · ffff256
1 parent 4ae5ae0
commit ffff256
Show file tree

Hide file tree

Showing 6 changed files with 107 additions and 103 deletions.
diff --git a/config/longrun_configs/longrun_moist_baroclinic_wave.yml b/config/longrun_configs/longrun_moist_baroclinic_wave.yml
@@ -4,6 +4,7 @@ dz_bottom: 30.0 # 300.0
 dt: "60secs" # 200
 approximate_linear_solve_iters: 3
 max_newton_iters_ode: 1
+rayleigh_sponge: true
 t_end: "10days"
 ode_algo: ARS343
 initial_condition: "MoistBaroclinicWave"

diff --git a/examples/Manifest-v1.11.toml b/examples/Manifest-v1.11.toml
@@ -359,7 +359,9 @@ weakdeps = ["CUDA", "MPI"]
 
 [[deps.ClimaCore]]
 deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "StaticArrays", "Statistics", "Unrolled"]
-git-tree-sha1 = "61d071cfe584d99a1400a3bf96a122c45c282121"
+git-tree-sha1 = "ce9de8f8c74e523a560b66b067c220438b287974"
+repo-rev = "main"
+repo-url = "https://github.com/CliMA/ClimaCore.jl.git"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
 version = "0.14.24"
 weakdeps = ["CUDA", "Krylov"]
@@ -418,8 +420,8 @@ version = "0.1.1"
 
 [[deps.ClimaTimeSteppers]]
 deps = ["ClimaComms", "Colors", "DataStructures", "DiffEqBase", "KernelAbstractions", "Krylov", "LinearAlgebra", "LinearOperators", "NVTX", "SciMLBase", "StaticArrays"]
-git-tree-sha1 = "f03e9f4316d380cdf851ec2c4c55efbfdb064439"
-repo-rev = "dy/move_T_imp"
+git-tree-sha1 = "facfe93769c3343bbe512fb61b19cfdf85804118"
+repo-rev = "main"
 repo-url = "https://github.com/CliMA/ClimaTimeSteppers.jl.git"
 uuid = "595c0a79-7f3d-439a-bc5a-b232dc3bde79"
 version = "0.8.1"

diff --git a/src/cache/precipitation_precomputed_quantities.jl b/src/cache/precipitation_precomputed_quantities.jl
@@ -43,19 +43,20 @@ function set_sedimentation_precomputed_quantities!(Y, p, t)
 
     (; ᶜwₗ, ᶜwᵢ) = p.precomputed
     cmc = CAP.microphysics_cloud_params(p.params)
+    FT = eltype(p.params)
 
     # compute the precipitation terminal velocity [m/s]
-    @. ᶜwₗ = CMNe.terminal_velocity(
-        cmc.liquid,
-        cmc.Ch2022.rain,
-        Y.c.ρ,
-        max(zero(Y.c.ρ), Y.c.ρq_liq / Y.c.ρ),
-    )
-    @. ᶜwᵢ = CMNe.terminal_velocity(
-        cmc.ice,
-        cmc.Ch2022.small_ice,
-        Y.c.ρ,
-        max(zero(Y.c.ρ), Y.c.ρq_ice / Y.c.ρ),
-    )
+    @. ᶜwₗ = FT(0) #CMNe.terminal_velocity(
+    #    cmc.liquid,
+    #    cmc.Ch2022.rain,
+    #    Y.c.ρ,
+    #    max(zero(Y.c.ρ), Y.c.ρq_liq / Y.c.ρ),
+    #)
+    @. ᶜwᵢ = FT(0) #CMNe.terminal_velocity(
+    #    cmc.ice,
+    #    cmc.Ch2022.small_ice,
+    #    Y.c.ρ,
+    #    max(zero(Y.c.ρ), Y.c.ρq_ice / Y.c.ρ),
+    #)
     return nothing
 end
diff --git a/src/cache/precomputed_quantities.jl b/src/cache/precomputed_quantities.jl
@@ -542,21 +542,21 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     @. ᶜwₕhₜ = Geometry.WVector(0)
     if moisture_model isa NonEquilMoistModel
         set_sedimentation_precomputed_quantities!(Y, p, t)
-        (; ᶜwₗ, ᶜwᵢ) = p.precomputed
-        @. ᶜwₜqₜ += Geometry.WVector(ᶜwₗ * Y.c.ρq_liq + ᶜwᵢ * Y.c.ρq_ice) / Y.c.ρ
-        @. ᶜwₕhₜ += Geometry.WVector(
-            ᶜwₗ * Y.c.ρq_liq * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwₗ) + Geometry.UVWVector(ᶜu))/2) +
-            ᶜwᵢ * Y.c.ρq_ice * (TD.internal_energy_ice(thermo_params, ᶜts)    + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwᵢ) + Geometry.UVWVector(ᶜu))/2)
-        ) / Y.c.ρ
+        #(; ᶜwₗ, ᶜwᵢ) = p.precomputed
+        #@. ᶜwₜqₜ += Geometry.WVector(ᶜwₗ * Y.c.ρq_liq + ᶜwᵢ * Y.c.ρq_ice) / Y.c.ρ
+        #@. ᶜwₕhₜ += Geometry.WVector(
+        #    ᶜwₗ * Y.c.ρq_liq * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwₗ) + Geometry.UVWVector(ᶜu))/2) +
+        #    ᶜwᵢ * Y.c.ρq_ice * (TD.internal_energy_ice(thermo_params, ᶜts)    + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwᵢ) + Geometry.UVWVector(ᶜu))/2)
+        #) / Y.c.ρ
     end
     if precip_model isa Microphysics1Moment
         set_precipitation_precomputed_quantities!(Y, p, t)
         (; ᶜwᵣ, ᶜwₛ) = p.precomputed
-        @. ᶜwₜqₜ += Geometry.WVector(ᶜwᵣ * Y.c.ρq_rai + ᶜwₛ * Y.c.ρq_sno) / Y.c.ρ
-        @. ᶜwₕhₜ += Geometry.WVector(
-            ᶜwᵣ * Y.c.ρq_rai * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwᵣ) + Geometry.UVWVector(ᶜu))/2) +
-            ᶜwₛ * Y.c.ρq_sno * (TD.internal_energy_ice(thermo_params, ᶜts)    + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwₛ) + Geometry.UVWVector(ᶜu))/2)
-        ) / Y.c.ρ
+        #@. ᶜwₜqₜ += Geometry.WVector(ᶜwᵣ * Y.c.ρq_rai + ᶜwₛ * Y.c.ρq_sno) / Y.c.ρ
+        #@. ᶜwₕhₜ += Geometry.WVector(
+        #    ᶜwᵣ * Y.c.ρq_rai * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwᵣ) + Geometry.UVWVector(ᶜu))/2) +
+        #    ᶜwₛ * Y.c.ρq_sno * (TD.internal_energy_ice(thermo_params, ᶜts)    + ᶜΦ - norm_sqr(Geometry.UVWVector(0, 0, ᶜwₛ) + Geometry.UVWVector(ᶜu))/2)
+        #) / Y.c.ρ
     end
 
     if turbconv_model isa PrognosticEDMFX

diff --git a/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl b/src/parameterized_tendencies/microphysics/microphysics_wrappers.jl
@@ -173,54 +173,54 @@ function compute_precipitation_sources!(
     @. Sqᵢᵖ -= Sᵖ
     @. Sqₛᵖ += Sᵖ
 
-    # accretion: q_liq + q_sno -> q_sno or q_rai
-    # sink of cloud water via accretion cloud water + snow
-    @. Sᵖ = min(
-        limit(qₗ(thp, ts, qᵣ), dt, 5),
-        CM1.accretion(mp.cl, mp.ps, mp.tv.snow, mp.ce, qₗ(thp, ts, qᵣ), qₛ, ρ),
-    )
-    # if T < T_freeze cloud droplets freeze to become snow
-    # else the snow melts and both cloud water and snow become rain
-    α(thp, ts) = TD.Parameters.cv_l(thp) / TD.latent_heat_fusion(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze)
-    @. Sᵖ_snow = ifelse(
-        Tₐ(thp, ts) < mp.ps.T_freeze,
-        Sᵖ,
-        FT(-1) * min(Sᵖ * α(thp, ts), limit(qₛ, dt, 5)),
-    )
-    @. Sqₛᵖ += Sᵖ_snow
-    @. Sqₗᵖ -= Sᵖ
-    @. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow)
-
-    # accretion: q_ice + q_rai -> q_sno
-    @. Sᵖ = min(
-        limit(qᵢ(thp, ts, qₛ), dt, 5),
-        CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₛ), qᵣ, ρ),
-    )
-    @. Sqᵢᵖ -= Sᵖ
-    @. Sqₛᵖ += Sᵖ
-    # sink of rain via accretion cloud ice - rain
-    @. Sᵖ = min(
-        limit(qᵣ, dt, 5),
-        CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₛ), qᵣ, ρ),
-    )
-    @. Sqᵣᵖ -= Sᵖ
-    @. Sqₛᵖ += Sᵖ
-
-    # accretion: q_rai + q_sno -> q_rai or q_sno
-    @. Sᵖ = ifelse(
-        Tₐ(thp, ts) < mp.ps.T_freeze,
-        min(
-            limit(qᵣ, dt, 5),
-            CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₛ, qᵣ, ρ),
-        ),
-        -min(
-            limit(qₛ, dt, 5),
-            CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qᵣ, qₛ, ρ),
-        ),
-    )
-    @. Sqₛᵖ += Sᵖ
-    @. Sqᵣᵖ -= Sᵖ
-    #! format: on
+    ## accretion: q_liq + q_sno -> q_sno or q_rai
+    ## sink of cloud water via accretion cloud water + snow
+    #@. Sᵖ = min(
+    #    limit(qₗ(thp, ts, qᵣ), dt, 5),
+    #    CM1.accretion(mp.cl, mp.ps, mp.tv.snow, mp.ce, qₗ(thp, ts, qᵣ), qₛ, ρ),
+    #)
+    ## if T < T_freeze cloud droplets freeze to become snow
+    ## else the snow melts and both cloud water and snow become rain
+    #α(thp, ts) = TD.Parameters.cv_l(thp) / TD.latent_heat_fusion(thp, ts) * (Tₐ(thp, ts) - mp.ps.T_freeze)
+    #@. Sᵖ_snow = ifelse(
+    #    Tₐ(thp, ts) < mp.ps.T_freeze,
+    #    Sᵖ,
+    #    FT(-1) * min(Sᵖ * α(thp, ts), limit(qₛ, dt, 5)),
+    #)
+    #@. Sqₛᵖ += Sᵖ_snow
+    #@. Sqₗᵖ -= Sᵖ
+    #@. Sqᵣᵖ += ifelse(Tₐ(thp, ts) < mp.ps.T_freeze, FT(0), Sᵖ - Sᵖ_snow)
+
+    ## accretion: q_ice + q_rai -> q_sno
+    #@. Sᵖ = min(
+    #    limit(qᵢ(thp, ts, qₛ), dt, 5),
+    #    CM1.accretion(mp.ci, mp.pr, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₛ), qᵣ, ρ),
+    #)
+    #@. Sqᵢᵖ -= Sᵖ
+    #@. Sqₛᵖ += Sᵖ
+    ## sink of rain via accretion cloud ice - rain
+    #@. Sᵖ = min(
+    #    limit(qᵣ, dt, 5),
+    #    CM1.accretion_rain_sink(mp.pr, mp.ci, mp.tv.rain, mp.ce, qᵢ(thp, ts, qₛ), qᵣ, ρ),
+    #)
+    #@. Sqᵣᵖ -= Sᵖ
+    #@. Sqₛᵖ += Sᵖ
+
+    ## accretion: q_rai + q_sno -> q_rai or q_sno
+    #@. Sᵖ = ifelse(
+    #    Tₐ(thp, ts) < mp.ps.T_freeze,
+    #    min(
+    #        limit(qᵣ, dt, 5),
+    #        CM1.accretion_snow_rain(mp.ps, mp.pr, mp.tv.rain, mp.tv.snow, mp.ce, qₛ, qᵣ, ρ),
+    #    ),
+    #    -min(
+    #        limit(qₛ, dt, 5),
+    #        CM1.accretion_snow_rain(mp.pr, mp.ps, mp.tv.snow, mp.tv.rain, mp.ce, qᵣ, qₛ, ρ),
+    #    ),
+    #)
+    #@. Sqₛᵖ += Sᵖ
+    #@. Sqᵣᵖ -= Sᵖ
+    ##! format: on
 end
 
 """
@@ -255,29 +255,29 @@ function compute_precipitation_sinks!(
     sps = (mp.ps, mp.tv.snow, mp.aps, thp)
     rps = (mp.pr, mp.tv.rain, mp.aps, thp)
 
-    #! format: off
-    # evaporation: q_rai -> q_vap
-    @. Sᵖ = -min(
-        limit(qᵣ, dt, 5),
-        -CM1.evaporation_sublimation(rps..., PP(thp, ts), qᵣ, ρ, Tₐ(thp, ts)),
-    )
-    @. Sqᵣᵖ += Sᵖ
-
-    # melting: q_sno -> q_rai
-    @. Sᵖ = min(
-        limit(qₛ, dt, 5),
-        CM1.snow_melt(sps..., qₛ, ρ, Tₐ(thp, ts)),
-    )
-    @. Sqᵣᵖ += Sᵖ
-    @. Sqₛᵖ -= Sᵖ
-
-    # deposition/sublimation: q_vap <-> q_sno
-    @. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₛ, ρ, Tₐ(thp, ts))
-    @. Sᵖ = ifelse(
-        Sᵖ > FT(0),
-        min(limit(qᵥ(thp, ts), dt, 5), Sᵖ),
-        -min(limit(qₛ, dt, 5), FT(-1) * Sᵖ),
-    )
-    @. Sqₛᵖ += Sᵖ
-    #! format: on
+    ##! format: off
+    ## evaporation: q_rai -> q_vap
+    #@. Sᵖ = -min(
+    #    limit(qᵣ, dt, 5),
+    #    -CM1.evaporation_sublimation(rps..., PP(thp, ts), qᵣ, ρ, Tₐ(thp, ts)),
+    #)
+    #@. Sqᵣᵖ += Sᵖ
+
+    ## melting: q_sno -> q_rai
+    #@. Sᵖ = min(
+    #    limit(qₛ, dt, 5),
+    #    CM1.snow_melt(sps..., qₛ, ρ, Tₐ(thp, ts)),
+    #)
+    #@. Sqᵣᵖ += Sᵖ
+    #@. Sqₛᵖ -= Sᵖ
+
+    ## deposition/sublimation: q_vap <-> q_sno
+    #@. Sᵖ = CM1.evaporation_sublimation(sps..., PP(thp, ts), qₛ, ρ, Tₐ(thp, ts))
+    #@. Sᵖ = ifelse(
+    #    Sᵖ > FT(0),
+    #    min(limit(qᵥ(thp, ts), dt, 5), Sᵖ),
+    #    -min(limit(qₛ, dt, 5), FT(-1) * Sᵖ),
+    #)
+    #@. Sqₛᵖ += Sᵖ
+    ##! format: on
 end
diff --git a/src/prognostic_equations/hyperdiffusion.jl b/src/prognostic_equations/hyperdiffusion.jl
@@ -244,11 +244,11 @@ NVTX.@annotate function apply_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     # TODO: Figure out why caching the duplicated tendencies in ᶜtemp_scalar
     # triggers allocations.
     for (ᶜρχₜ, ᶜ∇²χ, χ_name) in matching_subfields(Yₜ.c, ᶜ∇²specific_tracers)
-        #ν₄_scalar = ifelse(χ_name in (:q_rai, :q_sno), 0 * ν₄_scalar, ν₄_scalar)
+        ν₄_scalar = ifelse(χ_name in (:q_rai, :q_sno), 0 * ν₄_scalar, ν₄_scalar)
         @. ᶜρχₜ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
-        #if !(χ_name in (:q_rai, :q_sno))
-        @. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
-        #end
+        if !(χ_name in (:q_rai, :q_sno))
+            @. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
+        end
     end
     if turbconv_model isa PrognosticEDMFX
         (; ᶜ∇²q_totʲs) = p.hyperdiff