add convergence plot

experiments/standalone/Vegetation/timestep_test.jl

@@ -1,10 +1,47 @@
+# # Timestep test for standalone CanopyModel
+
+# This experiment runs the standalone canopy model for 100 days at a single point
+# and outputs useful plots related to the timestepping and stability of
+# the simulation using different timesteps. This was initially implemented to
+# test the behavior of stepping canopy temperature (T) implicitly.
+
+# The simulation is run multiple times: first, with a small reference timestep,
+# then with increasing timestep sizes. The goal here is to see which timesteps
+# are small enough to produce a stable simulation, and which timesteps are too
+# large to produce accurate results.
+
+# The runs with larger timesteps are compared with the reference timestep run
+# using multiple metrics, including:
+# - mean, 99%th, and 95th% error over the course of the simulation
+# - T at the end of the simulation (used to test convergence with each dt)
+# - T throughout the entire simulation
+
+# Note that the simulation is run for 100 days + 80 seconds. This somewhat
+# unusual length is necessary to test the convergence behavior of the solver.
+# Convergence behavior must be assessed when T is actively changing, rather
+# than when it is in equilibrium. T changes in the time period after a driver
+# update, which occurs every 3 hours (so exactly at 100 days). Running for
+# 100 days + 80 seconds allows us to look at results over a longer simulation,
+# as well as convergence behavior. See the discussion in
+# https://github.com/CliMA/ClimaLand.jl/pull/675 for more information.
+
+# Simulation Setup
+# Space: single point
+# Simulation duration: 100 days + 80 seconds (= 2400.02 hours)
+# Timestep: variable, ranging from 0.6s to 3600s
+# Timestepper: ARS111
+# Fixed number of iterations: 6
+# Jacobian update: Every Newton iteration
+# Atmospheric and radiation driver updates: Every 3 hours
+
 import SciMLBase
 import ClimaComms
 @static pkgversion(ClimaComms) >= v"0.6" && ClimaComms.@import_required_backends
 using CairoMakie
 using Statistics
 using Dates
 using Insolation
+using DelimitedFiles
 
 # Load CliMA Packages and ClimaLand Modules:
 
@@ -77,16 +114,12 @@ rt_params = TwoStreamParameters(
     λ_γ_PAR = FT(5e-7),
     λ_γ_NIR = FT(1.65e-6),
 )
-
 rt_model = TwoStreamModel{FT}(rt_params);
 
 cond_params = MedlynConductanceParameters(FT; g1 = FT(141.0))
-
 stomatal_model = MedlynConductanceModel{FT}(cond_params);
 
-
 photo_params = FarquharParameters(FT, Canopy.C3(); Vcmax25 = FT(5e-5))
-
 photosynthesis_model = FarquharModel{FT}(photo_params);
 
 AR_params = AutotrophicRespirationParameters(FT)
@@ -150,18 +183,13 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}(;
     radiation = radiation,
 );
 
-
-Y, p, coords = ClimaLand.initialize(canopy)
 exp_tendency! = make_exp_tendency(canopy)
 imp_tendency! = make_imp_tendency(canopy)
 jacobian! = make_jacobian(canopy);
-jac_kwargs =
-    (; jac_prototype = ClimaLand.ImplicitEquationJacobian(Y), Wfact = jacobian!);
-
+drivers = ClimaLand.get_drivers(canopy)
 
 ψ_leaf_0 = FT(-2e5 / 9800)
 ψ_stem_0 = FT(-1e5 / 9800)
-
 S_l_ini =
     inverse_water_retention_curve.(
         retention_model,
@@ -170,25 +198,16 @@ S_l_ini =
         S_s,
     )
 
-Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1])
-Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2])
-
-
 seconds_per_day = 3600 * 24.0
 t0 = 150seconds_per_day
-N_days = 100
-tf = t0 + N_days * seconds_per_day
-evaluate!(Y.canopy.energy.T, atmos.T, t0)
+N_days = 100.0
+# End simulation just after 100 days to see effect of driver update at 100 days
+tf = t0 + N_days * seconds_per_day + 80
+sim_time = round((tf - t0) / 3600, digits = 2) # simulation length in hours
 set_initial_cache! = make_set_initial_cache(canopy)
-set_initial_cache!(p, Y, t0);
-
-saveat = Array(t0:(3 * 3600):tf)
-updateat = Array(t0:(3600 * 3):tf)
-drivers = ClimaLand.get_drivers(canopy)
-updatefunc = ClimaLand.make_update_drivers(drivers)
-cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
 
 timestepper = CTS.ARS111();
+err = (FT == Float64) ? 1e-8 : 1e-4
 ode_algo = CTS.IMEXAlgorithm(
     timestepper,
     CTS.NewtonsMethod(
@@ -197,53 +216,112 @@ ode_algo = CTS.IMEXAlgorithm(
     ),
 );
 
-prob = SciMLBase.ODEProblem(
-    CTS.ClimaODEFunction(
-        T_exp! = exp_tendency!,
-        T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...),
-        dss! = ClimaLand.dss!,
-    ),
-    Y,
-    (t0, tf),
-    p,
-);
-
 ref_dt = 6.0
-ref_sol =
-    SciMLBase.solve(prob, ode_algo; dt = ref_dt, callback = cb, saveat = saveat);
-ref_T = [parent(ref_sol.u[k].canopy.energy.T)[1] for k in 1:length(ref_sol.t)]
+dts = [ref_dt, 12.0, 48.0, 225.0, 450.0, 900.0, 1800.0, 3600.0]
+
 mean_err = []
 p95_err = []
 p99_err = []
-dts = [12.0, 24.0, 48.0, 100.0, 225.0, 450.0, 900.0, 1800.0, 3600.0]
+sol_err = []
+T_states = []
+times = []
 for dt in dts
     @info dt
-    saveat = Array(t0:(3 * 3600):tf)
+
+    # Initialize model before each simulation
+    Y, p, coords = ClimaLand.initialize(canopy)
+    jac_kwargs = (;
+        jac_prototype = ClimaLand.ImplicitEquationJacobian(Y),
+        Wfact = jacobian!,
+    )
+
+    # Set initial conditions
+    Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1])
+    Y.canopy.hydraulics.ϑ_l.:2 .= augmented_liquid_fraction.(ν, S_l_ini[2])
     evaluate!(Y.canopy.energy.T, atmos.T, t0)
-    updateat = Array(t0:(3600 * 3):tf)
+    set_initial_cache!(p, Y, t0)
+
+    saveat = vcat(Array(t0:(3 * 3600):tf), tf)
+    updateat = vcat(Array(t0:(3 * 3600):tf), tf)
     updatefunc = ClimaLand.make_update_drivers(drivers)
     cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc)
 
+    prob = SciMLBase.ODEProblem(
+        CTS.ClimaODEFunction(
+            T_exp! = exp_tendency!,
+            T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...),
+            dss! = ClimaLand.dss!,
+        ),
+        Y,
+        (t0, tf),
+        p,
+    )
+
     @time sol =
         SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat)
-    T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]
-    ΔT = abs.(T .- ref_T)
-    push!(mean_err, mean(ΔT))
-    push!(p95_err, percentile(ΔT, 95))
-    push!(p99_err, percentile(ΔT, 99))
+    if dt == ref_dt
+        ref_T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]
+    else
+        T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)]
+
+        ΔT = abs.(T .- ref_T)
+        push!(mean_err, mean(ΔT))
+        push!(p95_err, percentile(ΔT, 95))
+        push!(p99_err, percentile(ΔT, 99))
+        push!(sol_err, T[end] - ref_T[end])
+        push!(T_states, T)
+        push!(times, sol.t)
+    end
 end
+
 savedir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Vegetation");
+
+# Create plot with statistics
+# Compare T state computed with small vs large dt
+mean_err = FT(mean(ΔT))
+p95_err = FT(percentile(ΔT, 95))
+p99_err = FT(percentile(ΔT, 99))
 fig = Figure()
 ax = Axis(
     fig[1, 1],
-    xlabel = "Time (minutes)",
+    xlabel = "Timestep (minutes)",
     ylabel = "Temperature (K)",
     xscale = log10,
     yscale = log10,
+    title = "Error in T over $(sim_time / 24) day sim, dts in [$(dts[1]), $(dts[end])]",
 )
 dts = dts ./ 60
 lines!(ax, dts, FT.(mean_err), label = "Mean Error")
 lines!(ax, dts, FT.(p95_err), label = "95th% Error")
 lines!(ax, dts, FT.(p99_err), label = "99th% Error")
 axislegend(ax)
 save(joinpath(savedir, "errors.png"), fig)
+
+# Create convergence plot
+fig2 = Figure()
+ax2 = Axis(
+    fig2[1, 1],
+    xlabel = "log(dt)",
+    ylabel = "log(|T[end] - T_ref[end]|)",
+    xscale = log10,
+    yscale = log10,
+    title = "Convergence of T; sim time = $(sim_time / 24) days, dts in [$(dts[1]), $(dts[end])]",
+)
+scatter!(ax2, dts, FT.(sol_err))
+lines!(ax2, dts, dts / 1e3)
+save(joinpath(savedir, "convergence.png"), fig2)
+
+# Plot T throughout full simulation for each run
+fig3 = Figure()
+ax3 = Axis(
+    fig3[1, 1],
+    xlabel = "time (hr)",
+    ylabel = "T (K)",
+    title = "T throughout simulation; length = $(sim_time / 24) days, dts in [$(dts[1]), $(dts[end])]",
+)
+times = times ./ 3600.0 # hours
+for i in 1:length(times)
+    lines!(ax3, times[i], T_states[i], label = "dt $(dts[i]) s")
+end
+axislegend(ax3, position = :rt)
+save(joinpath(savedir, "states.png"), fig3)

test/standalone/Vegetation/canopy_model.jl

@@ -236,8 +236,10 @@ import ClimaParams
             imp_tendency! = ClimaLand.make_imp_tendency(canopy)
             jacobian! = ClimaLand.make_jacobian(canopy)
             # set up jacobian info
-            jac_kwargs =
-                (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!)
+            jac_kwargs = (;
+                jac_prototype = ImplicitEquationJacobian(Y),
+                Wfact = jacobian!,
+            )
 
             t0 = FT(0.0)
             dY = similar(Y)
@@ -760,8 +762,10 @@ end
             imp_tendency! = ClimaLand.make_imp_tendency(canopy)
             jacobian! = ClimaLand.make_jacobian(canopy)
             # set up jacobian info
-            jac_kwargs =
-                (; jac_prototype = ImplicitEquationJacobian(Y), Wfact = jacobian!)
+            jac_kwargs = (;
+                jac_prototype = ImplicitEquationJacobian(Y),
+                Wfact = jacobian!,
+            )
 
             dY = similar(Y)
             exp_tendency!(dY, Y, p, t0)