canopy implicit timestepper

CliMA · Aug 16, 2024 · 8323271 · 8323271
1 parent 7c9bc04
commit 8323271
Show file tree

Hide file tree

Showing 20 changed files with 1,219 additions and 417 deletions.
diff --git a/.buildkite/Manifest.toml b/.buildkite/Manifest.toml
diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml
@@ -134,6 +134,15 @@ steps:
         env:
           CLIMACOMMS_DEVICE: "CUDA"
 
+      - label: "Canopy Implicit Stepping GPU"
+        command: "julia --color=yes --project=.buildkite experiments/standalone/Vegetation/timestep_test.jl"
+        artifact_paths: "experiments/standalone/Vegetation/errors.png"
+        agents:
+          slurm_ntasks: 1
+          slurm_gres: "gpu:1"
+        env:
+          CLIMACOMMS_DEVICE: "CUDA"
+
   - group: "ClimaLandSimulations"
     steps:
       - label: "Ozark figures Makie"

diff --git a/docs/Manifest.toml b/docs/Manifest.toml
diff --git a/experiments/benchmarks/land.jl b/experiments/benchmarks/land.jl
@@ -603,7 +603,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"
@@ -612,12 +612,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
@@ -629,7 +629,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
@@ -347,7 +347,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
@@ -356,7 +356,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
@@ -370,7 +378,7 @@ if drivers.GPP.status == absent
     plot_daily_avg(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         num_days,
         "μmol/m^2/s",
         savedir,
@@ -382,7 +390,7 @@ else
     plot_avg_comp(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         GPP_data,
         FT(DATA_DT),
         num_days,
@@ -397,7 +405,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,
@@ -410,7 +418,7 @@ else
     plot_avg_comp(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         SW_out_data,
         FT(DATA_DT),
         num_days,
@@ -425,7 +433,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,
@@ -438,7 +446,7 @@ else
     plot_avg_comp(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         LW_out_data,
         FT(DATA_DT),
         num_days,
@@ -460,7 +468,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) .*
@@ -469,7 +485,7 @@ else
     plot_avg_comp(
         "ET",
         ET_model,
-        dt * n,
+        dt_save,
         ET_data,
         FT(DATA_DT),
         num_days,
@@ -489,7 +505,7 @@ if drivers.H.status == absent
     plot_daily_avg(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -500,7 +516,7 @@ else
     plot_avg_comp(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         SHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
@@ -517,13 +533,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
@@ -291,7 +291,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
@@ -300,7 +300,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
@@ -314,7 +322,7 @@ if drivers.GPP.status == absent
     plot_daily_avg(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         num_days,
         "μmol/m^2/s",
         savedir,
@@ -326,7 +334,7 @@ else
     plot_avg_comp(
         "GPP",
         model_GPP,
-        dt * n,
+        dt_save,
         GPP_data,
         FT(DATA_DT),
         num_days,
@@ -341,7 +349,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,
@@ -354,7 +362,7 @@ else
     plot_avg_comp(
         "SW OUT",
         SW_out_model,
-        dt * n,
+        dt_save,
         SW_out_data,
         FT(DATA_DT),
         num_days,
@@ -369,7 +377,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,
@@ -382,7 +390,7 @@ else
     plot_avg_comp(
         "LW OUT",
         LW_out_model,
-        dt * n,
+        dt_save,
         LW_out_data,
         FT(DATA_DT),
         num_days,
@@ -404,7 +412,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) .*
@@ -413,7 +429,7 @@ else
     plot_avg_comp(
         "ET",
         ET_model,
-        dt * n,
+        dt_save,
         ET_data,
         FT(DATA_DT),
         num_days,
@@ -433,7 +449,7 @@ if drivers.H.status == absent
     plot_daily_avg(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         num_days,
         "w/m^2",
         savedir,
@@ -444,7 +460,7 @@ else
     plot_avg_comp(
         "SHF",
         SHF_model,
-        dt * n,
+        dt_save,
         SHF_data,
         FT(DATA_DT),
         N_days - N_spinup_days,
@@ -461,13 +477,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,