Add a GPU longrun with 1M microphysics (aquaplanet grid mean)

.buildkite/longruns_gpu/pipeline.yml

@@ -153,7 +153,7 @@ steps:
           slurm_time: 12:00:00
         env:
           JOB_NAME: "longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_diagedmf_0M"
-      
+
       - label: ":computer: aquaplanet equilmoist allsky radiation + diagnostic edmf + 0M microphysics"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
@@ -173,7 +173,7 @@ steps:
           slurm_time: 12:00:00
         env:
           JOB_NAME: "longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_0M_earth"
-      
+
       - label: ":computer: aquaplanet equilmoist clearsky radiation + 0M microphysics + earth topography (SLEVE)"
         command:
           - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
@@ -184,6 +184,16 @@ steps:
         env:
           JOB_NAME: "longrun_aquaplanet_rhoe_equil_55km_nz63_clearsky_0M_earth_sleve"
 
+      - label: ":umbrella: aquaplanet equilmoist clearsky radiation + 1M microphysics"
+        command:
+          - srun julia --project=examples examples/hybrid/driver.jl --config_file $CONFIG_PATH/$$JOB_NAME.yml
+        artifact_paths: "$$JOB_NAME/*"
+        agents:
+          slurm_gpus: 1
+          slurm_time: 6:00:00
+        env:
+          JOB_NAME: "longrun_aquaplanet_clearsky_1M"
+
   - group: "DYAMOND"
 
     steps:

config/longrun_configs/longrun_aquaplanet_clearsky_1M.yml

@@ -0,0 +1,18 @@
+dt: "120secs"
+t_end: "100days"
+h_elem: 16
+z_elem: 63
+dz_bottom: 30.0
+dz_top: 3000.0
+z_max: 55000.0
+rayleigh_sponge: true
+moist: "equil"
+precip_model: "1M"
+surface_setup: "DefaultMoninObukhov"
+vert_diff: "FriersonDiffusion"
+implicit_diffusion: true
+approximate_linear_solve_iters: 2
+rad: "clearsky"
+dt_rad: "6hours"
+job_id: "gpu_aquaplanet_clearsky_1M"
+toml: [toml/gpu_aquaplanet_clearsky_1M.toml]

docs/src/equations.md

@@ -385,13 +385,15 @@ It is assummed that some fraction ``\alpha`` of snow is melted during the proces
 ### Stability and positivity
 
 All source terms are individually limited such that they don't exceed the
-  available tracer specific humidity.
+  available tracer specific humidity divided by a coefficient ``a``.
 ```math
-\mathcal{S}_{x \rightarrow y} = min(\mathcal{S}_{x \rightarrow y}, \frac{q_{x}}{dt})
+\mathcal{S}_{x \rightarrow y} = min(\mathcal{S}_{x \rightarrow y}, \frac{q_{x}}{a \; dt})
 ```
 This will not ensure positivity because the sum of all source terms,
   combined with the advection tendency,
   could still drive the solution to negative numbers.
 It should however help mitigate some of the problems.
 The source terms functions treat negative specific humidities as zeros,
   so the simulations should be stable even with small negative numbers.
+
+We do not apply hyperdiffusion for precipitation tracers.

src/parameterized_tendencies/microphysics/precipitation.jl

@@ -190,7 +190,7 @@ function qₚ(ρqₚ::FT, ρ::FT) where {FT}
 end
 
 function limit(q::FT, dt::FT) where {FT}
-    return q / dt
+    return q / dt / FT(5)
 end
 
 function compute_precipitation_cache!(Y, p, colidx, ::Microphysics1Moment, _)

src/prognostic_equations/hyperdiffusion.jl

@@ -236,7 +236,7 @@ NVTX.@annotate function apply_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     h_length_scale = Spaces.node_horizontal_length_scale(h_space) # mean nodal distance
     ν₄_scalar = ν₄_scalar_coeff * h_length_scale^3
     n = n_mass_flux_subdomains(turbconv_model)
-
+    FT = eltype(Y)
     (; ᶜ∇²specific_tracers) = p.hyperdiff
 
     # TODO: Since we are not applying the limiter to density (or area-weighted
@@ -246,7 +246,7 @@ NVTX.@annotate function apply_tracer_hyperdiffusion_tendency!(Yₜ, Y, p, t)
     # triggers allocations.
     for (ᶜρχₜ, ᶜ∇²χ, χ_name) in matching_subfields(Yₜ.c, ᶜ∇²specific_tracers)
         ν₄_scalar =
-            ifelse(χ_name in (:q_rai, :q_sno), 0.1 * ν₄_scalar, ν₄_scalar)
+            ifelse(χ_name in (:q_rai, :q_sno), FT(0) * ν₄_scalar, ν₄_scalar)
         @. ᶜρχₜ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))
         if !(χ_name in (:q_rai, :q_sno))
             @. Yₜ.c.ρ -= ν₄_scalar * wdivₕ(Y.c.ρ * gradₕ(ᶜ∇²χ))