canopy implicit timestepper

CliMA · Sep 25, 2024 · e6c817d · e6c817d
1 parent 735590e
commit e6c817d
Show file tree

Hide file tree

Showing 21 changed files with 1,478 additions and 118 deletions.
diff --git a/.buildkite/pipeline.yml b/.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"

diff --git a/docs/src/APIs/canopy/Photosynthesis.md b/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

diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
@@ -606,7 +606,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"
@@ -615,12 +615,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
@@ -632,7 +632,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[]

diff --git a/experiments/integrated/fluxnet/US-Ha1/US-Ha1_simulation.jl b/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)
diff --git a/experiments/integrated/fluxnet/US-MOz/US-MOz_simulation.jl b/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)
diff --git a/experiments/integrated/fluxnet/US-NR1/US-NR1_simulation.jl b/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)
diff --git a/experiments/integrated/fluxnet/US-Var/US-Var_simulation.jl b/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)
diff --git a/experiments/integrated/fluxnet/fluxnet_simulation.jl b/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),
     ),
 );
diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl
@@ -348,7 +348,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
@@ -357,7 +357,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
@@ -371,7 +379,7 @@ if drivers.GPP.status == absent
     plot_daily_avg(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         num_days,
         "μmol/m^2/s",
         savedir,
@@ -383,7 +391,7 @@ else
     plot_avg_comp(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         GPP_data,
         FT(DATA_DT),
         num_days,
@@ -398,7 +406,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,
@@ -411,7 +419,7 @@ else
     plot_avg_comp(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         SW_out_data,
         FT(DATA_DT),
         num_days,
@@ -426,7 +434,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,
@@ -439,7 +447,7 @@ else
     plot_avg_comp(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         LW_out_data,
         FT(DATA_DT),
         num_days,
@@ -463,7 +471,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) .*
@@ -472,7 +488,7 @@ else
     plot_avg_comp(
         "ET",
         ET_model,
-        dt * n,
+        dt_save,
         ET_data,
         FT(DATA_DT),
         num_days,
@@ -492,7 +508,7 @@ if drivers.H.status == absent
     plot_daily_avg(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -503,7 +519,7 @@ else
     plot_avg_comp(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         SHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
@@ -520,13 +536,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,

diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl
@@ -292,7 +292,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
@@ -301,7 +301,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
@@ -315,7 +323,7 @@ if drivers.GPP.status == absent
     plot_daily_avg(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         num_days,
         "μmol/m^2/s",
         savedir,
@@ -327,7 +335,7 @@ else
     plot_avg_comp(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         GPP_data,
         FT(DATA_DT),
         num_days,
@@ -342,7 +350,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,
@@ -355,7 +363,7 @@ else
     plot_avg_comp(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         SW_out_data,
         FT(DATA_DT),
         num_days,
@@ -370,7 +378,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,
@@ -383,7 +391,7 @@ else
     plot_avg_comp(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         LW_out_data,
         FT(DATA_DT),
         num_days,
@@ -407,7 +415,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) .*
@@ -416,7 +432,7 @@ else
     plot_avg_comp(
         "ET",
         ET_model,
-        dt * n,
+        dt_save,
         ET_data,
         FT(DATA_DT),
         num_days,
@@ -436,7 +452,7 @@ if drivers.H.status == absent
     plot_daily_avg(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -447,7 +463,7 @@ else
     plot_avg_comp(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         SHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
@@ -464,13 +480,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,