Merge pull request #675 from CliMA/kd/canopy_temp_implicit

Add implicit tendency for canopy temperature

.buildkite/Manifest.toml

@@ -371,9 +371,9 @@ weakdeps = ["CUDA", "MPI"]
 
 [[deps.ClimaCore]]
 deps = ["Adapt", "BandedMatrices", "BlockArrays", "ClimaComms", "CubedSphere", "DataStructures", "DocStringExtensions", "ForwardDiff", "GaussQuadrature", "GilbertCurves", "HDF5", "InteractiveUtils", "IntervalSets", "KrylovKit", "LinearAlgebra", "MultiBroadcastFusion", "NVTX", "PkgVersion", "RecursiveArrayTools", "RootSolvers", "SparseArrays", "StaticArrays", "Statistics", "Unrolled"]
-git-tree-sha1 = "13c2f4e58c78fa54a22705d15e039c99462112ed"
+git-tree-sha1 = "6db06e8234dc67c021d636936803300ea710e0c6"
 uuid = "d414da3d-4745-48bb-8d80-42e94e092884"
-version = "0.14.17"
+version = "0.14.18"
 weakdeps = ["CUDA", "Krylov"]
 
     [deps.ClimaCore.extensions]

.buildkite/pipeline.yml

@@ -127,6 +127,14 @@ steps:
         agents:
           slurm_mem: 16G
 
+      - label: "Canopy Implicit Stepping CPU"
+        command: "julia --color=yes --project=.buildkite experiments/standalone/Vegetation/timestep_test.jl"
+        artifact_paths: "experiments/standalone/Vegetation/timestep_test/*"
+
+      - label: "SoilCanopy Implicit Stepping CPU"
+        command: "julia --color=yes --project=.buildkite experiments/integrated/performance/integrated_timestep_test.jl"
+        artifact_paths: "experiments/integrated/performance/integrated_timestep_test/*"
+
   - group: "Experiments on GPU"
     steps:
       - label: "Richards Runoff GPU"

NEWS.md

@@ -6,6 +6,7 @@ main
 - Add a longrun simulation for the bucket PR[#807](https://github.com/CliMA/ClimaLand.jl/pull/807)
 - Add ground heat flux to snow-soil model PR[#796](https://github.com/CliMA/ClimaLand.jl/pull/796)
 - Add snow-soil model PR [#779](https://github.com/CliMA/ClimaLand.jl/pull/779)
+- Step canopy temperature implicitly. PR [#675](https://github.com/CliMA/ClimaLand.jl/pull/675)
 
 v0.15.0
 --------

docs/src/APIs/canopy/Photosynthesis.md

@@ -20,8 +20,6 @@ ClimaLand.Canopy.intercellular_co2
 ClimaLand.Canopy.co2_compensation
 ClimaLand.Canopy.rubisco_assimilation
 ClimaLand.Canopy.light_assimilation
-ClimaLand.Canopy.C3
-ClimaLand.Canopy.C4
 ClimaLand.Canopy.max_electron_transport
 ClimaLand.Canopy.electron_transport
 ClimaLand.Canopy.net_photosynthesis

experiments/benchmarks/land.jl

@@ -608,7 +608,7 @@ function setup_simulation(; greet = false)
     Δt = 900.0
     nelements = (101, 15)
     if greet
-        @info "Run: Global RichardsModel"
+        @info "Run: Global Land Model"
         @info "Resolution: $nelements"
         @info "Timestep: $Δt s"
         @info "Duration: $(tf - t0) s"
@@ -617,12 +617,12 @@ function setup_simulation(; greet = false)
     prob, cb = setup_prob(t0, tf, Δt; nelements)
 
     # Define timestepper and ODE algorithm
-    stepper = CTS.ARS343()
+    stepper = CTS.ARS111()
     ode_algo = CTS.IMEXAlgorithm(
         stepper,
         CTS.NewtonsMethod(
-            max_iters = 1,
-            update_j = CTS.UpdateEvery(CTS.NewTimeStep),
+            max_iters = 6,
+            update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
         ),
     )
     return prob, ode_algo, Δt, cb
@@ -634,7 +634,7 @@ SciMLBase.solve(prob, ode_algo; dt = Δt, callback = cb)
 
 @info "Starting profiling"
 # Stop when we profile for MAX_PROFILING_TIME_SECONDS or MAX_PROFILING_SAMPLES
-MAX_PROFILING_TIME_SECONDS = 1000
+MAX_PROFILING_TIME_SECONDS = 500
 MAX_PROFILING_SAMPLES = 100
 time_now = time()
 timings_s = Float64[]

experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl

@@ -20,7 +20,4 @@ h_stem = FT(14) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(40)
-
-# Number of timesteps between saving output:
-n = 45
+dt = Float64(450)

experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl

@@ -20,7 +20,4 @@ h_leaf = FT(9.5) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(450)
-
-# Number of timesteps between saving output:
-n = 4
+dt = Float64(900)

experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl

@@ -20,7 +20,4 @@ h_stem = FT(7.5) # m
 t0 = Float64(120 * 3600 * 24)# start mid year to avoid snow
 
 # Time step size:
-dt = Float64(40)
-
-# Number of timesteps between saving output:
-n = 45
+dt = Float64(450)

experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl

@@ -20,7 +20,4 @@ h_stem = FT(0) # m
 t0 = Float64(21 * 3600 * 24)# start day 21 of the year
 
 # Time step size:
-dt = Float64(20)
-
-# Number of timesteps between saving output:
-n = 45
+dt = Float64(900)

experiments/integrated/fluxnet/fluxnet_simulation.jl

@@ -6,14 +6,15 @@ N_spinup_days = 30
 N_days = N_spinup_days + 30
 tf = Float64(t0 + 3600 * 24 * N_days)
 t_spinup = Float64(t0 + N_spinup_days * 3600 * 24)
-saveat = Array(t_spinup:(n * dt):tf)
+dt_save = 3600.0
+saveat = Array(t_spinup:dt_save:tf)
 
 # Set up timestepper
-timestepper = CTS.ARS343();
+timestepper = CTS.ARS111();
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
-        max_iters = 1,
+        max_iters = 6,
         update_j = CTS.UpdateEvery(CTS.NewNewtonIteration),
     ),
 );

experiments/integrated/fluxnet/ozark_pft.jl

@@ -344,7 +344,7 @@ num_days = N_days - N_spinup_days
 
 # Time series of model and data outputs
 data_times = [0:DATA_DT:(num_days * S_PER_DAY);]
-model_times = [0:(n * dt):(num_days * S_PER_DAY);]
+model_times = [0:dt_save:(num_days * S_PER_DAY);]
 
 # Plot model diurnal cycles without data comparisons
 # Autotrophic Respiration
@@ -353,7 +353,15 @@ AR =
         parent(sv.saveval[k].canopy.autotrophic_respiration.Ra)[1] for
         k in 1:length(sv.saveval)
     ] .* 1e6
-plot_daily_avg("AutoResp", AR, dt * n, num_days, "μmol/m^2/s", savedir, "Model")
+plot_daily_avg(
+    "AutoResp",
+    AR,
+    dt_save,
+    num_days,
+    "μmol/m^2/s",
+    savedir,
+    "Model",
+)
 
 # Plot all comparisons of model diurnal cycles to data diurnal cycles
 # GPP
@@ -367,7 +375,7 @@ if drivers.GPP.status == absent
     plot_daily_avg(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         num_days,
         "μmol/m^2/s",
         savedir,
@@ -379,7 +387,7 @@ else
     plot_avg_comp(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         GPP_data,
         FT(DATA_DT),
         num_days,
@@ -394,7 +402,7 @@ if drivers.SW_OUT.status == absent
     plot_daily_avg(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -407,7 +415,7 @@ else
     plot_avg_comp(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         SW_out_data,
         FT(DATA_DT),
         num_days,
@@ -422,7 +430,7 @@ if drivers.LW_OUT.status == absent
     plot_daily_avg(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -435,7 +443,7 @@ else
     plot_avg_comp(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         LW_out_data,
         FT(DATA_DT),
         num_days,
@@ -459,7 +467,15 @@ E =
     ] .* (1e3 * 24 * 3600)
 ET_model = T .+ E
 if drivers.LE.status == absent
-    plot_daily_avg("ET", ET_model, dt * n, num_days, "mm/day", savedir, "Model")
+    plot_daily_avg(
+        "ET",
+        ET_model,
+        dt_save,
+        num_days,
+        "mm/day",
+        savedir,
+        "Model",
+    )
 else
     measured_T =
         drivers.LE.values ./ (LP.LH_v0(earth_param_set) * 1000) .*
@@ -468,7 +484,7 @@ else
     plot_avg_comp(
         "ET",
         ET_model,
-        dt * n,
+        dt_save,
         ET_data,
         FT(DATA_DT),
         num_days,
@@ -488,7 +504,7 @@ if drivers.H.status == absent
     plot_daily_avg(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -499,7 +515,7 @@ else
     plot_avg_comp(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         SHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
@@ -516,13 +532,13 @@ LHF_canopy =
     [parent(sv.saveval[k].canopy.energy.lhf)[1] for k in 1:length(sol.t)]
 LHF_model = LHF_soil + LHF_canopy
 if drivers.LE.status == absent
-    plot_daily_avg("LHF", LHF_model, dt * n, num_days, "w/m^2", savedir)
+    plot_daily_avg("LHF", LHF_model, dt_save, num_days, "w/m^2", savedir)
 else
     LHF_data = drivers.LE.values[Int64(t_spinup ÷ DATA_DT):Int64(tf ÷ DATA_DT)]
     plot_avg_comp(
         "LHF",
         LHF_model,
-        dt * n,
+        dt_save,
         LHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,