diff --git a/docs/Manifest.toml b/docs/Manifest.toml index d03178cc8b..e9f11fbec5 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1,6 +1,6 @@ # This file is machine-generated - editing it directly is not advised -julia_version = "1.10.9" +julia_version = "1.10.10" manifest_format = "2.0" project_hash = "0e313427333375504d76ddc31cd582bee13918ac" @@ -2186,7 +2186,7 @@ version = "3.2.4+0" [[deps.OpenLibm_jll]] deps = ["Artifacts", "Libdl"] uuid = "05823500-19ac-5b8b-9628-191a04bc5112" -version = "0.8.1+4" +version = "0.8.5+0" [[deps.OpenMPI_jll]] deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "MPIPreferences", "TOML", "Zlib_jll"] diff --git a/docs/src/tutorials/integrated/soil_canopy_tutorial.jl b/docs/src/tutorials/integrated/soil_canopy_tutorial.jl index be1e59d523..f6cd4e755b 100644 --- a/docs/src/tutorials/integrated/soil_canopy_tutorial.jl +++ b/docs/src/tutorials/integrated/soil_canopy_tutorial.jl @@ -50,6 +50,7 @@ using ClimaLand.Soil using ClimaLand.Soil.Biogeochemistry using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand import ClimaLand.Parameters as LP using DelimitedFiles @@ -65,7 +66,11 @@ earth_param_set = LP.LandParameters(FT); site_ID = "US-MOz" # Timezone (offset from UTC in hrs) time_offset = 7 -start_date = DateTime(2010) + Hour(time_offset) +start_date = DateTime(2010) + Hour(time_offset) + Day(150) # start mid year +# Specify the time range and dt value over which to perform the simulation. +N_days = 100 +end_date = start_date + Day(N_days) +dt = Float64(30) # Site latitude and longitude lat = FT(38.7441) # degree @@ -74,11 +79,6 @@ long = FT(-92.2000) # degree # Height of the sensor at the site atmos_h = FT(32) # m -# Specify the time range and dt value over which to perform the simulation. -t0 = Float64(150 * 3600 * 24)# start mid year -N_days = 100 -tf = t0 + Float64(3600 * 24 * N_days) -dt = Float64(30) # - We will be using prescribed atmospheric and radiative drivers from the # US-MOz tower, which we read in here. We are using prescribed @@ -186,8 +186,8 @@ soilco2 = Soil.Biogeochemistry.SoilCO2Model{FT}(domain, drivers); # Read in prescribed LAI at the site from global MODIS data surface_space = domain.space.surface; modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(; - start_date = start_date + Second(t0), - end_date = start_date + Second(t0) + Second(tf), + start_date = start_date, + end_date = end_date, ); LAI = ClimaLand.prescribed_lai_modis( modis_lai_ncdata_path, @@ -273,50 +273,40 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}( # Now we can combine the soil and canopy models into a single combined model. land = SoilCanopyModel{FT}(soilco2, soil, canopy); -# Now we can initialize the state vectors and model coordinates, and initialize -# the explicit/implicit tendencies as usual. The Richard's equation time -# stepping is done implicitly, while the canopy model may be explicitly stepped, -# so we use an IMEX (implicit-explicit) scheme for the combined model. - -Y, p, coords = initialize(land); -exp_tendency! = make_exp_tendency(land); -imp_tendency! = make_imp_tendency(land); -jacobian! = make_jacobian(land); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); - # We need to provide initial conditions for the model: -Y.soil.ϑ_l = FT(0.4) -Y.soil.θ_i = FT(0.0) -T_0 = FT(288.7) -ρc_s = - volumetric_heat_capacity.( - Y.soil.ϑ_l, - Y.soil.θ_i, - land.soil.parameters.ρc_ds, - earth_param_set, - ) -Y.soil.ρe_int = - volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T_0, earth_param_set) - -Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air - -canopy_retention_model = canopy.hydraulics.parameters.retention_model -ψ_stem_0 = FT(-1e5 / 9800) -ψ_leaf_0 = FT(-2e5 / 9800) - -S_l_ini = - inverse_water_retention_curve.( - canopy_retention_model, - [ψ_stem_0, ψ_leaf_0], - plant_ν, - plant_S_s, - ) -for i in 1:2 - Y.canopy.hydraulics.ϑ_l.:($i) .= - augmented_liquid_fraction.(plant_ν, S_l_ini[i]) +function set_ic!(Y, p, t0, model) + Y.soil.ϑ_l = FT(0.4) + Y.soil.θ_i = FT(0.0) + T_0 = FT(288.7) + ρc_s = + volumetric_heat_capacity.( + Y.soil.ϑ_l, + Y.soil.θ_i, + land.soil.parameters.ρc_ds, + earth_param_set, + ) + Y.soil.ρe_int = + volumetric_internal_energy.(Y.soil.θ_i, ρc_s, T_0, earth_param_set) + + Y.soilco2.C .= FT(0.000412) # set to atmospheric co2, mol co2 per mol air + + canopy_retention_model = canopy.hydraulics.parameters.retention_model + ψ_stem_0 = FT(-1e5 / 9800) + ψ_leaf_0 = FT(-2e5 / 9800) + + S_l_ini = + inverse_water_retention_curve.( + canopy_retention_model, + [ψ_stem_0, ψ_leaf_0], + plant_ν, + plant_S_s, + ) + for i in 1:2 + Y.canopy.hydraulics.ϑ_l.:($i) .= + augmented_liquid_fraction.(plant_ν, S_l_ini[i]) + end + evaluate!(Y.canopy.energy.T, atmos.T, t0) end -evaluate!(Y.canopy.energy.T, atmos.T, t0) # Choose the timestepper and solver needed for the problem. timestepper = CTS.ARS343() @@ -328,47 +318,37 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Now set the initial values for the cache variables for the combined soil and plant model. -set_initial_cache! = make_set_initial_cache(land) -set_initial_cache!(p, Y, t0); -# Set the callbacks, which govern +# Set the callbacks update times, which govern # how often we save output, and how often we update # the forcing data ("drivers") n = 120 -saveat = Array(t0:(n * dt):tf) +saveat = Array(start_date:Second(n * dt):end_date); sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) -model_drivers = ClimaLand.get_drivers(land) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -updateat = Array(t0:1800:tf) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); - -# Carry out the simulation -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - callback = cb, - adaptive = false, - saveat = saveat, +updateat = Array(start_date:Second(1800):end_date); + +# Create the LandSimulation object, which will also create and initialize the state vectors, +# the cache, the driver callbacks, and set the initial conditions. +simulation = LandSimulation( + start_date, + end_date, + dt, + land; + set_ic! = set_ic!, + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), ); +sol = solve!(simulation); + # # Plotting # Now that we have both a soil and canopy model incorporated together, we will @@ -376,7 +356,7 @@ sol = SciMLBase.solve( # each of these models. As before, we may plot the GPP of the system as well # as transpiration showing fluxes in the canopy. -daily = sol.t ./ 3600 ./ 24 +daily = FT.(sol.t) ./ 3600 ./ 24 model_GPP = [ parent(sv.saveval[k].canopy.photosynthesis.GPP)[1] for k in 1:length(sv.saveval) diff --git a/docs/src/tutorials/standalone/Bucket/bucket_tutorial.jl b/docs/src/tutorials/standalone/Bucket/bucket_tutorial.jl index c9dbb2b98a..36ba962109 100644 --- a/docs/src/tutorials/standalone/Bucket/bucket_tutorial.jl +++ b/docs/src/tutorials/standalone/Bucket/bucket_tutorial.jl @@ -144,6 +144,7 @@ import ClimaParams as CP using ClimaLand.Bucket: BucketModel, BucketModelParameters, PrescribedBaregroundAlbedo using ClimaLand.Domains: coordinates, Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand: initialize, make_update_aux, @@ -203,10 +204,6 @@ z_0b = FT(1e-3); ρc_soil = FT(2e6); # Snow melt timescale τc = FT(3600); -# Simulation start date, end time, and timestep -t0 = 0.0; -tf = 7 * 86400; -Δt = 3600.0; bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc); @@ -215,6 +212,10 @@ bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc); # time, the date of the start of the simulation. In this tutorial we will # consider this January 1, 2005. start_date = DateTime(2005); +# Simulation start date, end time, and timestep + +end_date = start_date + Second(7 * 86400); +Δt = 3600.0; # To drive the system in standalone mode, # the user must provide @@ -230,9 +231,9 @@ start_date = DateTime(2005); # Precipitation: precip = (t) -> 0; -snow_precip = (t) -> -5e-7 * (t > 3 * 86400) * (t < 4 * 86400); +snow_precip = (t) -> -5e-7 * (FT(t) > 3 * 86400) * (FT(t) < 4 * 86400); # Diurnal temperature variations: -T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2); +T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * FT(t) / 86400 - π / 2); # Constant otherwise: u_atmos = (t) -> 3.0; q_atmos = (t) -> 0.005; @@ -256,8 +257,8 @@ bucket_atmos = PrescribedAtmosphere( # peak at local noon, and a prescribed downwelling LW radiative # flux, assuming the air temperature is on average 275 degrees # K with a diurnal amplitude of 5 degrees K: -SW_d = (t) -> @. max(1361 * sin(2π * t / 86400 - π / 2)); -LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2))^4; +SW_d = (t) -> @. max(1361 * sin(2π * FT(t) / 86400 - π / 2)); +LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * FT(t) / 86400 - π / 2))^4; bucket_rad = PrescribedRadiativeFluxes( FT, TimeVaryingInput(SW_d), @@ -281,89 +282,88 @@ model = BucketModel( # type and rely on multiple dispatch to obtain the atmospheric and radiative # quantitites from the coupler. -# Like all ClimaLand models, we set up the state vector using `initialize`: -Y, p, coords = initialize(model); - -# We can inspect the prognostic and auxiliary variables of the model: -ClimaLand.prognostic_vars(model) -Y.bucket |> propertynames -# The auxiliary variables in this case are the surface temperature, the turbulent fluxes, the -# net radiation, and the surface specific humidity. -ClimaLand.auxiliary_vars(model) -p.bucket |> propertynames -# Next is to set initial conditions. -Y.bucket.T .= FT(270); -Y.bucket.W .= FT(0.05); -Y.bucket.Ws .= FT(0.0); -Y.bucket.σS .= FT(0.08); -# We also set the initial values of the cache here: -set_initial_cache! = make_set_initial_cache(model); -set_initial_cache!(p, Y, t0); - -# Then to create the entire right hand side (tendency) function for the system -# of ordinary differential equations: -exp_tendency! = make_exp_tendency(model); +# Next is to create a function to set initial conditions. +function set_ic!(Y, p, t0, model) + Y.bucket.T .= FT(270) + Y.bucket.W .= FT(0.05) + Y.bucket.Ws .= FT(0.0) + Y.bucket.σS .= FT(0.08) +end # Now we choose our timestepping algorithm. timestepper = CTS.RK4() ode_algo = CTS.ExplicitAlgorithm(timestepper) -# Then we can set up the simulation and solve it: -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!), - Y, - (t0, tf), - p, -); + # We need a callback to get and store the auxiliary fields, as they # are not stored by default. We also need a callback to update the # drivers (atmos and radiation) -saveat = collect(t0:Δt:tf); +saveat = collect(start_date:Second(Δt):end_date); saved_values = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ); saving_cb = ClimaLand.NonInterpSavingCallback(saved_values, saveat); -updateat = copy(saveat) -model_drivers = ClimaLand.get_drivers(model) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb) +updateat = deepcopy(saveat); -sol = SciMLBase.solve(prob, ode_algo; dt = Δt, saveat = saveat, callback = cb); +# Create the LandSimulation object, which will also create and initialize the state vectors, +# the cache, the driver callbacks, and set the initial conditions. +simulation = LandSimulation( + start_date, + end_date, + Δt, + model; + set_ic! = set_ic!, + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), +); +Y = simulation._integrator.u; +p = simulation._integrator.p; +# We can inspect the prognostic and auxiliary variables of the model: +ClimaLand.prognostic_vars(model) +Y.bucket |> propertynames +# The auxiliary variables in this case are the surface temperature, the turbulent fluxes, the +# net radiation, and the surface specific humidity. +ClimaLand.auxiliary_vars(model) +p.bucket |> propertynames +sol = solve!(simulation); # Extracting the solution from what is returned by the ODE.jl commands # is a bit clunky right now, but we are working on hiding some of this. # `parent` extracts the underlying data from the ClimaCore.Fields.Field object # and we loop over the solution `sol` because of how the data is stored # within solutions returned by ODE.jl - indexed by timestep. -W = [parent(sol.u[k].bucket.W)[1] for k in 1:length(sol.t)]; -Ws = [parent(sol.u[k].bucket.Ws)[1] for k in 1:length(sol.t)]; -σS = [parent(sol.u[k].bucket.σS)[1] for k in 1:length(sol.t)]; +times = FT.(sol.t) +W = [parent(sol.u[k].bucket.W)[1] for k in 1:length(times)]; +Ws = [parent(sol.u[k].bucket.Ws)[1] for k in 1:length(times)]; +σS = [parent(sol.u[k].bucket.σS)[1] for k in 1:length(times)]; T_sfc = - [parent(saved_values.saveval[k].bucket.T_sfc)[1] for k in 1:length(sol.t)]; + [parent(saved_values.saveval[k].bucket.T_sfc)[1] for k in 1:length(times)]; evaporation = [ parent(saved_values.saveval[k].bucket.turbulent_fluxes.vapor_flux)[1] - for k in 1:length(sol.t) + for k in 1:length(times) ]; -R_n = [parent(saved_values.saveval[k].bucket.R_n)[1] for k in 1:length(sol.t)]; +R_n = [parent(saved_values.saveval[k].bucket.R_n)[1] for k in 1:length(times)]; # The turbulent energy flux is the sum of latent and sensible heat fluxes. LHF = [ parent(saved_values.saveval[k].bucket.turbulent_fluxes.lhf)[1] for - k in 1:length(sol.t) + k in 1:length(times) ]; SHF = [ parent(saved_values.saveval[k].bucket.turbulent_fluxes.shf)[1] for - k in 1:length(sol.t) + k in 1:length(times) ]; turbulent_energy_flux = SHF .+ LHF plot( - sol.t ./ 86400, + times ./ 86400, W, label = "", xlabel = "time (days)", @@ -374,7 +374,7 @@ savefig("w.png") # ![](w.png) plot( - sol.t ./ 86400, + times ./ 86400, σS, label = "", xlabel = "time (days)", @@ -385,20 +385,20 @@ savefig("swe.png") # ![](swe.png) plot( - sol.t ./ 86400, - snow_precip.(sol.t), + times ./ 86400, + snow_precip.(times), label = "Net precipitation", xlabel = "time (days)", ylabel = "Flux (m/s)", title = "Surface water fluxes", legend = :bottomright, ) -plot!(sol.t ./ 86400, evaporation, label = "Sublimation/Evaporation") +plot!(times ./ 86400, evaporation, label = "Sublimation/Evaporation") savefig("water_f.png") # ![](water_f.png) plot( - sol.t ./ 86400, + times ./ 86400, T_sfc, title = "Surface Temperatures", label = "Ground temperature", @@ -406,11 +406,11 @@ plot( ylabel = "T_sfc (K)", legend = :bottomright, ) -plot!(sol.t ./ 86400, T_atmos.(sol.t), label = "Atmospheric Temperature") +plot!(times ./ 86400, T_atmos.(times), label = "Atmospheric Temperature") savefig("t.png") # ![](t.png) plot( - sol.t ./ 86400, + times ./ 86400, R_n, label = "Net radiative flux", xlabel = "time (days)", @@ -418,8 +418,8 @@ plot( title = "Surface energy fluxes", legend = :bottomright, ) -plot!(sol.t ./ 86400, turbulent_energy_flux, label = "Turbulent fluxes") -plot!(sol.t ./ 86400, R_n .+ turbulent_energy_flux, label = "Net flux") +plot!(times ./ 86400, turbulent_energy_flux, label = "Turbulent fluxes") +plot!(times ./ 86400, R_n .+ turbulent_energy_flux, label = "Net flux") savefig("energy_f.png") # ![](energy_f.png) diff --git a/docs/src/tutorials/standalone/Canopy/canopy_tutorial.jl b/docs/src/tutorials/standalone/Canopy/canopy_tutorial.jl index c99b174e76..d1882ec6c8 100644 --- a/docs/src/tutorials/standalone/Canopy/canopy_tutorial.jl +++ b/docs/src/tutorials/standalone/Canopy/canopy_tutorial.jl @@ -50,6 +50,7 @@ using ClimaLand using ClimaLand.Domains: Point using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand import ClimaLand.Parameters as LP using DelimitedFiles @@ -66,6 +67,13 @@ earth_param_set = LP.LandParameters(FT); time_offset = 7 start_date = DateTime(2010) + Hour(time_offset) +# Select a time range to perform time stepping over, and a dt. As usual, +# the timestep depends on the problem you are solving, the accuracy of the +# solution required, and the timestepping algorithm you are using. +N_days = 364 +end_date = start_date + Day(N_days) +dt = 225.0 + # Site latitude and longitude lat = FT(38.7441) # degree long = FT(-92.2000) # degree @@ -92,13 +100,7 @@ site_ID = "US-MOz"; # would change. domain = Point(; z_sfc = FT(0.0), longlat = (long, lat)); -# Select a time range to perform time stepping over, and a dt. As usual, -# the timestep depends on the problem you are solving, the accuracy of the -# solution required, and the timestepping algorithm you are using. -t0 = 0.0 -N_days = 364 -tf = t0 + 3600 * 24 * N_days -dt = 225.0; + # We will be using prescribed atmospheric and radiative drivers from the # US-MOz tower, which we read in here. We are using prescribed @@ -131,10 +133,8 @@ forcing = (; atmos, radiation, ground); # Now we read in time-varying LAI from a global MODIS dataset. surface_space = domain.space.surface; -modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths( - start_date = start_date + Second(t0), - end_date = start_date + Second(t0) + Second(tf), -) +modis_lai_ncdata_path = + ClimaLand.Artifacts.modis_lai_multiyear_paths(; start_date, end_date) LAI = ClimaLand.prescribed_lai_modis( modis_lai_ncdata_path, surface_space, @@ -208,17 +208,6 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}( photosynthesis, ); -# Initialize the state vectors and obtain the model coordinates, then get the -# explicit time stepping tendency that updates auxiliary and prognostic -# variables that are stepped explicitly. - -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.FieldMatrixWithSolver(Y), Wfact = jacobian!); - # Provide initial conditions for the canopy hydraulics model (; retention_model, ν, S_s) = canopy.hydraulics.parameters ψ_stem_0 = FT(-1e5 / 9800) @@ -231,35 +220,28 @@ S_l_ini = ν, S_s, ) -for i in 1:2 - Y.canopy.hydraulics.ϑ_l.:($i) .= augmented_liquid_fraction.(ν, S_l_ini[i]) +function set_ic!(Y, p, t0, model) + for i in 1:2 + Y.canopy.hydraulics.ϑ_l.:($i) .= + augmented_liquid_fraction.(ν, S_l_ini[i]) + end + evaluate!(Y.canopy.energy.T, atmos.T, t0) end -evaluate!(Y.canopy.energy.T, atmos.T, t0) - -# Initialize the cache variables for the canopy using the initial -# conditions and initial time. - -set_initial_cache! = make_set_initial_cache(canopy) -set_initial_cache!(p, Y, t0); # Allocate the struct which stores the saved auxiliary state # and create the callback which saves it at each element in saveat. n = 16 -saveat = Array(t0:(n * dt):tf) +saveat = Array(start_date:Second(n * dt):end_date); sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); # Create the callback function which updates the forcing variables, # or drivers. -updateat = Array(t0:1800:tf) -model_drivers = ClimaLand.get_drivers(canopy) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); +updateat = Array(start_date:Second(1800):end_date); # Select a timestepping algorithm and setup the ODE problem. @@ -272,28 +254,27 @@ 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, +# Create the LandSimulation object, which will also create and initialize the state vectors, +# the cache, the driver callbacks, and set the initial conditions. +simulation = LandSimulation( + start_date, + end_date, + dt, + canopy; + set_ic! = set_ic!, + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), ); -# Now, we can solve the problem and store the model data in the saveat array, -# using [`SciMLBase.jl`](https://github.com/SciML/SciMLBase.jl) and -# [`ClimaTimeSteppers.jl`](https://github.com/CliMA/ClimaTimeSteppers.jl). - -sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); +sol = solve!(simulation); # # Create some plots # We can now plot the data produced in the simulation. For example, GPP: -daily = sol.t ./ 3600 ./ 24 +daily = FT.(sol.t) ./ 3600 ./ 24 model_GPP = [ parent(sv.saveval[k].canopy.photosynthesis.GPP)[1] for k in 1:length(sv.saveval) diff --git a/docs/src/tutorials/standalone/Soil/evaporation.jl b/docs/src/tutorials/standalone/Soil/evaporation.jl index c27fd1ef9c..2b9aeb1d5c 100644 --- a/docs/src/tutorials/standalone/Soil/evaporation.jl +++ b/docs/src/tutorials/standalone/Soil/evaporation.jl @@ -17,6 +17,7 @@ using DelimitedFiles: readdlm using ClimaLand using ClimaLand.Domains: Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand.Soil import ClimaLand import ClimaLand.Parameters as LP @@ -39,7 +40,11 @@ thermo_params = LP.thermodynamic_parameters(earth_param_set); # Our assumption is that in the lab experiment there was no radiative heating # or cooling of the soil. -start_date = DateTime(2005) # required argument, but not used in this case +# Timestepping: +start_date = DateTime(2005) +end_date = start_date + Day(13) +dt = Float64(900.0) + SW_d = (t) -> 0 LW_d = (t) -> 301.15^4 * 5.67e-8 radiation = PrescribedRadiativeFluxes( @@ -132,7 +137,9 @@ function hydrostatic_equilibrium(z, z_interface, params) return ν * (1 + (α * (z - z_interface))^n)^(-m) end end -function init_soil!(Y, z, params) +function set_ic!(Y, p, t0, model) + params = model.parameters + z = model.domain.fields.z FT = eltype(Y.soil.ϑ_l) Y.soil.ϑ_l .= hydrostatic_equilibrium.(z, FT(-0.001), params) Y.soil.θ_i .= 0 @@ -162,25 +169,7 @@ soil = Soil.EnergyHydrology{FT}(; sources = (), ) -Y, p, cds = initialize(soil); - -init_soil!(Y, z, soil.parameters); - # Timestepping: -t0 = Float64(0) -tf = Float64(24 * 3600 * 13) -dt = Float64(900.0) - -# We also set the initial conditions of the cache here: -set_initial_cache! = make_set_initial_cache(soil) -set_initial_cache!(p, Y, t0); - -# Define the tendency functions -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!) timestepper = CTS.ARS111(); ode_algo = CTS.IMEXAlgorithm( @@ -191,32 +180,29 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Define the problem and callbacks: -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -saveat = Array(t0:3600.0:tf) +saveat = Array(start_date:Second(3600.0):end_date); sv_hr = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv_hr, saveat) updateat = deepcopy(saveat) -model_drivers = ClimaLand.get_drivers(soil) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); + +simulation = LandSimulation( + start_date, + end_date, + dt, + soil; + set_ic! = set_ic!, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + updateat = updateat, + diagnostics = (), +); # Solve -sol_hr = - SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); +sol_hr = solve!(simulation); # Repeat at lower resolution zmax = FT(0) @@ -233,26 +219,7 @@ soil = Soil.EnergyHydrology{FT}(; sources = (), ) -Y, p, cds = initialize(soil); - -init_soil!(Y, z, soil.parameters); - # Timestepping: -t0 = Float64(0) -tf = Float64(24 * 3600 * 13) -dt = Float64(3600.0) - -# We also set the initial conditions of the cache here: -set_initial_cache! = make_set_initial_cache(soil) -set_initial_cache!(p, Y, t0); - -# Define the tendency functions -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); - timestepper = CTS.ARS111(); ode_algo = CTS.IMEXAlgorithm( timestepper, @@ -262,32 +229,28 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Define the problem and callbacks: -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -saveat = Array(t0:3600.0:tf) +saveat = Array(start_date:Second(3600.0):end_date); sv_lr = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv_lr, saveat) updateat = deepcopy(saveat) -model_drivers = ClimaLand.get_drivers(soil) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); +simulation = LandSimulation( + start_date, + end_date, + dt, + soil; + set_ic! = set_ic!, + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), +); # Solve -sol_lr = - SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); +sol_lr = solve!(simulation); # Figures @@ -314,7 +277,7 @@ ax = Axis( xgridvisible = false, ygridvisible = false, ) -CairoMakie.xlims!(minimum(data_dates), maximum(sol_lr.t ./ 3600 ./ 24)) +CairoMakie.xlims!(minimum(data_dates), maximum(float.(sol_lr.t) ./ 3600 ./ 24)) CairoMakie.lines!( ax, FT.(data_dates), @@ -325,7 +288,7 @@ CairoMakie.lines!( ) CairoMakie.lines!( ax, - sol_lr.t ./ 3600 ./ 24, + FT.(sol_lr.t) ./ 3600 ./ 24, evap_lr .* (1000 * 3600 * 24), label = "Model, 7 elements", color = :blue, @@ -333,7 +296,7 @@ CairoMakie.lines!( ) CairoMakie.lines!( ax, - sol_hr.t ./ 3600 ./ 24, + FT.(sol_hr.t) ./ 3600 ./ 24, evap_hr .* (1000 * 3600 * 24), label = "Model, 28 elements", color = :blue, diff --git a/docs/src/tutorials/standalone/Soil/evaporation_gilat_loess.jl b/docs/src/tutorials/standalone/Soil/evaporation_gilat_loess.jl index 285dfa1ece..f679ca9072 100644 --- a/docs/src/tutorials/standalone/Soil/evaporation_gilat_loess.jl +++ b/docs/src/tutorials/standalone/Soil/evaporation_gilat_loess.jl @@ -29,6 +29,7 @@ using ClimaLand.Domains: Column using ClimaLand.Soil import ClimaLand import ClimaLand.Parameters as LP +import ClimaLand.Simulations: LandSimulation, solve! import SurfaceFluxes.Parameters as SFP FT = Float64; @@ -68,6 +69,8 @@ params = ClimaLand.Soil.EnergyHydrologyParameters( ); start_date = DateTime(2005) +end_date = start_date + Day(15) +dt = Float64(900.0) SW_d = (t) -> 0 LW_d = (t) -> 294.15^4 * 5.67e-8 radiation = PrescribedRadiativeFluxes( @@ -115,9 +118,6 @@ no_flux_boundary_fluxes = (; bottom = WaterHeatBC(; water = zero_water_flux, heat = zero_heat_flux), ); -t0 = Float64(0) -tf = Float64(24 * 3600 * 15) -dt = Float64(900.0) Δz = 0.01 zmax = FT(0) zmin = FT(-1.6) @@ -132,12 +132,12 @@ soil = Soil.EnergyHydrology{FT}(; sources = (), ); # Initial conditions - -Y, p, cds = initialize(soil) function estimated_ic(z) 0.34 / (1 + exp(-(z + 0.165) / 0.005)) + 0.05 end -function init_soil!(Y, z, params) +function set_ic!(Y, p, t0, model) + params = model.parameters + z = model.domain.fields.z FT = eltype(Y.soil.ϑ_l) Y.soil.ϑ_l .= estimated_ic.(z) Y.soil.θ_i .= 0 @@ -156,18 +156,6 @@ function init_soil!(Y, z, params) params.earth_param_set, ) end - -init_soil!(Y, z, soil.parameters) -set_initial_cache! = make_set_initial_cache(soil) -set_initial_cache!(p, Y, t0); - -# Timestepping: -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); - timestepper = CTS.ARS111(); ode_algo = CTS.IMEXAlgorithm( timestepper, @@ -177,27 +165,27 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Problem definition and callbacks -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -saveat = Array(t0:3600.0:tf) +# Saving callback +saveat = Array(start_date:Hour(1):end_date); sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) cb = SciMLBase.CallbackSet(saving_cb); - -sol_no_evap = - SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); +simulation = LandSimulation( + start_date, + end_date, + dt, + soil; + set_ic! = set_ic!, + updateat = [], + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), +); +sol_no_evap = solve!(simulation); # Repeat with evaporation and drainage # This requires different initial conditions @@ -213,46 +201,29 @@ soil = Soil.EnergyHydrology{FT}(; boundary_conditions = evap_boundary_fluxes, sources = (), ) -Y, p, cds = initialize(soil) -init_soil!(Y, z, soil.parameters) -set_initial_cache! = make_set_initial_cache(soil) -set_initial_cache!(p, Y, t0) -soil_exp_tendency! = make_exp_tendency(soil) -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); + timestepper = CTS.ARS111(); -ode_algo = CTS.IMEXAlgorithm( - timestepper, - CTS.NewtonsMethod( - max_iters = 1, - update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), - ), -); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -saveat = Array(t0:3600.0:tf) +saveat = Array(start_date:Hour(1):end_date); sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) updateat = deepcopy(saveat) -model_drivers = ClimaLand.get_drivers(soil) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb) -sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat) + +simulation = LandSimulation( + start_date, + end_date, + dt, + soil; + set_ic! = set_ic!, + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), +); +sol = solve!(simulation); evap = [ parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for k in 1:length(sol.t) @@ -279,46 +250,26 @@ soil = Soil.EnergyHydrology{FT}(; boundary_conditions = no_drainage_boundary_fluxes, sources = (), ) -Y, p, cds = initialize(soil) -init_soil!(Y, z_no_evap, soil.parameters) -set_initial_cache! = make_set_initial_cache(soil) -set_initial_cache!(p, Y, t0) -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -timestepper = CTS.ARS111(); -ode_algo = CTS.IMEXAlgorithm( - timestepper, - CTS.NewtonsMethod( - max_iters = 1, - update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), - ), -); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -saveat = Array(t0:3600.0:tf) +saveat = Array(start_date:Hour(1):end_date); sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) updateat = deepcopy(saveat) -model_drivers = ClimaLand.get_drivers(soil) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb) -sol_no_drainage = - SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat) +simulation = LandSimulation( + start_date, + end_date, + dt, + soil; + set_ic! = set_ic!, + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), +); +sol_no_drainage = solve!(simulation); evap_no_drainage = [ parent(sv.saveval[k].soil.turbulent_fluxes.vapor_flux_liq)[1] for k in 1:length(sol.t) @@ -331,14 +282,14 @@ fig = Figure(size = (800, 400)) ax = Axis(fig[1, 1], xlabel = "Day", ylabel = "Evaporation rate (mm/d)") CairoMakie.lines!( ax, - sol.t ./ 3600 ./ 24, + FT.(sol.t) ./ 3600 ./ 24, evap .* (1000 * 3600 * 24), label = "With drainage", color = :red, ) CairoMakie.lines!( ax, - sol_no_drainage.t ./ 3600 ./ 24, + FT.(sol_no_drainage.t) ./ 3600 ./ 24, evap_no_drainage .* (1000 * 3600 * 24), label = "No drainage", color = :blue, @@ -348,13 +299,13 @@ CairoMakie.axislegend(ax) ax2 = Axis(fig[1, 2], xlabel = "Day", ylabel = "Cumulative evaporation (mm)") CairoMakie.lines!( ax2, - sol.t ./ 3600 ./ 24, + FT.(sol.t) ./ 3600 ./ 24, cumsum(evap) .* (1000 * 3600), color = :red, ) CairoMakie.lines!( ax2, - sol_no_drainage.t ./ 3600 ./ 24, + FT.(sol_no_drainage.t) ./ 3600 ./ 24, cumsum(evap_no_drainage) .* (1000 * 3600), color = :blue, ) diff --git a/docs/src/tutorials/standalone/Soil/freezing_front.jl b/docs/src/tutorials/standalone/Soil/freezing_front.jl index 2fb2ddfbc8..a60e373964 100644 --- a/docs/src/tutorials/standalone/Soil/freezing_front.jl +++ b/docs/src/tutorials/standalone/Soil/freezing_front.jl @@ -82,6 +82,7 @@ import ClimaParams as CP using ClimaLand using ClimaLand.Domains: Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand.Soil import ClimaLand @@ -169,15 +170,11 @@ soil = Soil.EnergyHydrology{FT}(; # # Running a simulation -# Once we have the model, we can -# [`initialize`](@ref ClimaLand.initialize) -# the state vectors and obtain the coordinates -Y, p, coords = initialize(soil); +# We can specify the initial condition: -# After which, we can specify the initial condition -# function, and initialze the variables: - -function init_soil!(Ysoil, z, params) +function set_ic!(Ysoil, p, t0, model) + params = model.parameters + z = model.domain.fields.z ν = params.ν FT = eltype(Ysoil.soil.ϑ_l) Ysoil.soil.ϑ_l .= FT(0.33) @@ -198,23 +195,10 @@ function init_soil!(Ysoil, z, params) ) end -init_soil!(Y, coords.subsurface.z, soil.parameters); - # We choose the initial and final simulation times: t0 = Float64(0) tf = Float64(60 * 60 * 50); -# We set the cache values corresponding to the initial conditions -# of the state Y: -set_initial_cache! = make_set_initial_cache(soil); -set_initial_cache!(p, Y, t0); -# Create the tendency function, and choose a timestep, and timestepper: -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); - dt = Float64(100) timestepper = CTS.ARS111(); @@ -226,20 +210,20 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Problem definition and callbacks -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, +simulation = LandSimulation( + t0, + tf, + dt, + soil; + set_ic! = set_ic!, + solver_kwargs = (; save_everystep = true), + timestepper = ode_algo, + user_callbacks = (), + diagnostics = (), ); +sol = solve!(simulation); -# Now we can solve the problem. -sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = collect(0:3600:tf)); +# Now we can create and solve the simulation. # # Comparison to data # This data was obtained by us from the figures of Hansson et al. (2004), but was originally obtained @@ -289,7 +273,7 @@ ax3 = Axis( limits!(ax3, 0.2, 0.5, -0.2, 0.0) -z = parent(coords.subsurface.z)[:]; +z = parent(soil.domain.fields.z)[:]; scatter!(ax1, vwc[mask_12h], -depth[mask_12h], label = "", color = "orange") lines!( diff --git a/docs/src/tutorials/standalone/Soil/layered_soil.jl b/docs/src/tutorials/standalone/Soil/layered_soil.jl index c76c7357fc..6f4c486723 100644 --- a/docs/src/tutorials/standalone/Soil/layered_soil.jl +++ b/docs/src/tutorials/standalone/Soil/layered_soil.jl @@ -17,6 +17,7 @@ import ClimaParams as CP using DelimitedFiles: readdlm using ClimaLand +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand.Domains: Column using ClimaLand.Soil import ClimaLand @@ -97,15 +98,10 @@ soil = Soil.RichardsModel{FT}(; sources = (), ); -# Initial the state vectors, and set initial conditions -Y, p, cds = initialize(soil); - -# Initial conditions -Y.soil.ϑ_l .= 0.0353; # read from Figure 4 of Huang et al. - -# We also set the initial conditions of the auxiliary state here: -set_initial_cache! = make_set_initial_cache(soil) -set_initial_cache!(p, Y, t0); +# create function to initialize the state vectors +function set_ic!(Y, p, t0, model) + Y.soil.ϑ_l .= 0.0353 # read from Figure 4 of Huang et al. +end # Timestepping: stepper = CTS.ARS111() @@ -121,24 +117,19 @@ ode_algo = CTS.IMEXAlgorithm( convergence_checker = conv_checker, ), ) -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil) -jacobian! = make_jacobian(soil) - -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!) -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -) saveat = [0.0, 8.0, 16.0, 24.0, 32.0, 40.0, 60.0] .* 60 # chosen to compare with data in plots in paper -sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = saveat); +simulation = LandSimulation( + t0, + tf, + dt, + soil; + set_ic! = set_ic!, + solver_kwargs = (; saveat = saveat), + timestepper = ode_algo, + user_callbacks = (), + diagnostics = (), +); +sol = solve!(simulation); z = parent(ClimaCore.Fields.coordinate_field(soil_domain.space.subsurface).z) ϑ_l = [parent(sol.u[k].soil.ϑ_l) for k in 1:length(sol.t)] diff --git a/docs/src/tutorials/standalone/Soil/phase_change_analytic.jl b/docs/src/tutorials/standalone/Soil/phase_change_analytic.jl index 988f1b03f0..b9c395a483 100644 --- a/docs/src/tutorials/standalone/Soil/phase_change_analytic.jl +++ b/docs/src/tutorials/standalone/Soil/phase_change_analytic.jl @@ -57,6 +57,7 @@ import ClimaParams as CP using ClimaLand using ClimaLand.Domains: Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand.Soil import ClimaLand @@ -127,16 +128,13 @@ soil = Soil.EnergyHydrology{FT}(; # # Running a simulation -# Once we have the model, we can -# [`initialize`](@ref ClimaLand.initialize) -# the state vectors and obtain the coordinates -Y, p, coords = initialize(soil); -# After which, we can specify the initial condition -# function, and initialze the variables. We chose these to match +# Specify the initial condition +# function We chose these to match # the initial conditions of the Stefan problem: -function init_soil!(Ysoil, z, params) +function set_ic!(Ysoil, p, t0, model) + params = model.parameters ν = params.ν FT = eltype(Ysoil.soil.ϑ_l) Ysoil.soil.ϑ_l .= FT(0.33) @@ -157,23 +155,10 @@ function init_soil!(Ysoil, z, params) ) end -init_soil!(Y, coords.subsurface.z, soil.parameters); - # We choose the initial and final simulation times: t0 = Float64(0) tf = Float64(60 * 60 * 24 * 20); -# We set the cache values corresponding to the initial conditions -# of the state Y: -set_initial_cache! = make_set_initial_cache(soil); -set_initial_cache!(p, Y, t0); -# Create the tendency function, and choose a timestep, and timestepper: -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); - dt = Float64(100) timestepper = CTS.ARS111(); @@ -185,31 +170,28 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Problem definition and callbacks -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -saveat = Array(t0:3600.0:tf) +# Callbacks +saveat = Array(t0:3600.0:tf); sv = (; t = Array{Float64}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); -# Now we can solve the problem. -sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - saveat = collect(0:3600:tf), - callback = saving_cb, +# Now we can create and solve the simulation +simulation = LandSimulation( + t0, + tf, + dt, + soil; + set_ic! = set_ic!, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), ); + +# Solve +sol = solve!(simulation); sol_T = parent(sv.saveval[end].soil.T)[:] fig = Figure(size = (500, 500), fontsize = 24) @@ -224,7 +206,7 @@ ax1 = Axis( ) limits!(ax1, 262, 276, -3, 0.0) -z = parent(coords.subsurface.z)[:]; +z = parent(soil.domain.fields.z)[:]; lines!(ax1, sol_T, z, label = "Model", color = :blue, linewidth = 3) @@ -314,7 +296,7 @@ function implicit(ζ) end ζ = find_zero(implicit, (0.25, 0.27), Bisection()) depth = Array(0:0.01:3) -t = sol.t[end] +t = FT(sol.t[end]) zf = 2.0 * ζ * sqrt(d1 * t) analytic_unfrozen_profile(depth, zf) = erfc(depth / (zf / ζ / (d1 / d2)^0.5)) / erfc(ζ * (d1 / d2)^0.5) diff --git a/docs/src/tutorials/standalone/Soil/richards_equation.jl b/docs/src/tutorials/standalone/Soil/richards_equation.jl index 7ec2862cd7..e68c8003db 100644 --- a/docs/src/tutorials/standalone/Soil/richards_equation.jl +++ b/docs/src/tutorials/standalone/Soil/richards_equation.jl @@ -63,6 +63,7 @@ import ClimaTimeSteppers as CTS using ClimaLand using ClimaLand.Domains: Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand.Soil import ClimaLand @@ -135,37 +136,15 @@ soil = Soil.RichardsModel{FT}(; sources = sources, ); -# Here we create the explicit and implicit tendencies, which update prognostic -# variable components that are stepped explicitly and implicitly, respectively. -# We also create the function which is used to update our Jacobian. -exp_tendency! = make_exp_tendency(soil); -imp_tendency! = ClimaLand.make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); - -# # Set up the simulation -# We can now initialize the prognostic and auxiliary variable vectors, and take -# a peek at what those variables are: -Y, p, coords = initialize(soil); -Y.soil |> propertynames - -p.soil |> propertynames -coords |> propertynames - -# Note that the variables are nested into `Y` and `p` in a hierarchical way. -# Since we have the vectors (composed of [ClimaCore Fields](https://clima.github.io/ClimaCore.jl/dev/api/#ClimaCore.Fields.Field) handy, we can now set them to the desired initial -# conditions. -Y.soil.ϑ_l .= FT(0.494); +function set_ic!(Y, p, t0, model) + Y.soil.ϑ_l .= FT(0.494) +end # We choose the initial and final simulation times: t0 = Float64(0) tf = Float64(60 * 60 * 24 * 36); -# We set the cache values corresponding to the initial conditions -# of the state Y: -set_initial_cache! = make_set_initial_cache(soil); -set_initial_cache!(p, Y, t0); - # Next, we turn to timestepping. # As usual, your timestep depends on the problem you are solving, the accuracy # of the solution required, and the timestepping algorithm you are using. @@ -184,24 +163,31 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Here we set up the information used for our Jacobian. -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); - -# And then we can solve the system of equations, using -# [SciMLBase.jl](https://github.com/SciML/SciMLBase.jl) and -# [ClimaTimeSteppers.jl](https://github.com/CliMA/ClimaTimeSteppers.jl). -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, +simulation = LandSimulation( + t0, + tf, + dt, + soil; + set_ic! = set_ic!, + solver_kwargs = (; save_everystep = true), + timestepper = ode_algo, + user_callbacks = (), + diagnostics = (), ); -sol = SciMLBase.solve(prob, ode_algo; dt = dt, adaptive = false); + +# We can now take +# a peek at what the initialized variables are: +Y = simulation._integrator.u +p = simulation._integrator.p +Y.soil |> propertynames + +p.soil |> propertynames + + +# Note that the variables are nested into `Y` and `p` in a hierarchical way. +# Since we have the vectors (composed of [ClimaCore Fields](https://clima.github.io/ClimaCore.jl/dev/api/#ClimaCore.Fields.Field) handy, we can now set them to the desired initial +# conditions. +sol = solve!(simulation); # # Create some plots @@ -219,9 +205,9 @@ sol = SciMLBase.solve(prob, ode_algo; dt = dt, adaptive = false); # requiring that the integrated water content at the end matches that # at the beginning (yielding an interface location of `z≈-0.56m`). -t = sol.t ./ (60 * 60 * 24); +t = FT.(sol.t) ./ (60 * 60 * 24); ϑ_l = [parent(sol.u[k].soil.ϑ_l) for k in 1:length(t)] -z = parent(coords.subsurface.z) +z = parent(soil.domain.fields.z) plot( ϑ_l[1], z, diff --git a/docs/src/tutorials/standalone/Soil/soil_energy_hydrology.jl b/docs/src/tutorials/standalone/Soil/soil_energy_hydrology.jl index a289dca176..a542cd1268 100644 --- a/docs/src/tutorials/standalone/Soil/soil_energy_hydrology.jl +++ b/docs/src/tutorials/standalone/Soil/soil_energy_hydrology.jl @@ -88,6 +88,7 @@ import ClimaParams as CP import ClimaTimeSteppers as CTS using ClimaLand using ClimaLand.Domains: Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand.Soil import ClimaLand @@ -180,27 +181,11 @@ soil = Soil.EnergyHydrology{FT}(; sources = sources, ); -# Here we create the explicit and implicit tendencies, which update prognostic -# variable components that are stepped explicitly and implicitly, respectively. -# We also create the function which is used to update our Jacobian. -exp_tendency! = make_exp_tendency(soil); -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -# # Set up the simulation -# We can now initialize the prognostic and auxiliary variable vectors, and take -# a peek at what those variables are: -Y, p, coords = initialize(soil); -Y.soil |> propertynames -p.soil |> propertynames - -coords |> propertynames - -# Note that the variables are nested into `Y` and `p` in a hierarchical way. -# Since we have the vectors handy, we can now set them to the desired initial -# conditions. -function init_soil!(Y, z, params) +function set_ic!(Y, p, t0, model) + params = model.parameters + z = model.domain.fields.z ν = params.ν θ_r = params.θ_r FT = eltype(Y.soil.ϑ_l) @@ -237,17 +222,10 @@ function init_soil!(Y, z, params) ) end -init_soil!(Y, coords.subsurface.z, soil.parameters); - # We choose the initial and final simulation times: t0 = Float64(0) tf = Float64(60 * 60 * 72); -# We set the cache values corresponding to the initial conditions -# of the state Y: -set_initial_cache! = make_set_initial_cache(soil); -set_initial_cache!(p, Y, t0); - # We use [ClimaTimesteppers.jl](https://github.com/CliMA/ClimaTimesteppers.jl) for carrying out the time integration. # Choose a timestepper and set up the ODE problem: @@ -261,38 +239,35 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); - -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); - - # By default, it # only returns Y and t at each time we request output (`saveat`, below). We use # a callback in order to also get the auxiliary vector `p` back: -saveat = collect(t0:FT(30000):tf) +saveat = collect(t0:FT(30000):tf); saved_values = (; t = Array{Float64}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ); cb = ClimaLand.NonInterpSavingCallback(saved_values, saveat); +simulation = LandSimulation( + t0, + tf, + dt, + soil; + set_ic! = set_ic!, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (cb,), + diagnostics = (), +); + # Now we can solve the problem. -sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = saveat, callback = cb); +sol = solve!(simulation); # Extract output -z = parent(coords.subsurface.z) -t = parent(sol.t) +z = parent(soil.domain.fields.z) +t = parent(FT.(sol.t)) ϑ_l = [parent(sol.u[k].soil.ϑ_l) for k in 1:length(t)] T = [parent(saved_values.saveval[k].soil.T) for k in 1:length(t)]; # Let's look at the initial and final times: diff --git a/docs/src/tutorials/standalone/Soil/sublimation.jl b/docs/src/tutorials/standalone/Soil/sublimation.jl index b3467bbe7f..208fb32c79 100644 --- a/docs/src/tutorials/standalone/Soil/sublimation.jl +++ b/docs/src/tutorials/standalone/Soil/sublimation.jl @@ -17,6 +17,7 @@ using DelimitedFiles: readdlm using ClimaLand using ClimaLand.Domains: Column using ClimaLand.Soil +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand import ClimaLand.Parameters as LP import SurfaceFluxes.Parameters as SFP @@ -26,6 +27,8 @@ earth_param_set = LP.LandParameters(FT) thermo_params = LP.thermodynamic_parameters(earth_param_set); start_date = DateTime(2005) +end_date = start_date + Day(4) +dt = Float64(10) SW_d = (t) -> 0 LW_d = (t) -> 270.0^4 * 5.67e-8 radiation = PrescribedRadiativeFluxes( @@ -123,9 +126,7 @@ soil = Soil.EnergyHydrology{FT}(; sources = sources, ) -Y, p, cds = initialize(soil); - -# Set initial conditions +# functions to Set initial conditions function hydrostatic_equilibrium(z, z_interface, params) (; ν, S_s, hydrology_cm) = params (; α, n, m) = hydrology_cm @@ -135,8 +136,9 @@ function hydrostatic_equilibrium(z, z_interface, params) return ν * (1 + (α * (z - z_interface))^n)^(-m) end end -function init_soil!(Y, z, params) - FT = eltype(Y.soil.ϑ_l) +function set_ic!(Y, p, t0, model) + params = model.parameters + z = model.domain.fields.z Y.soil.ϑ_l .= hydrostatic_equilibrium.(z, FT(-0.1), params) Y.soil.θ_i .= 0 T = FT(275.0) @@ -149,24 +151,8 @@ function init_soil!(Y, z, params) Y.soil.ρe_int = Soil.volumetric_internal_energy.(FT(0), ρc_s, T, params.earth_param_set) end -init_soil!(Y, z, soil.parameters); # Timestepping: -t0 = Float64(0) -tf = Float64(24 * 3600 * 4) -dt = Float64(10) - -# We also set the initial conditions of the cache here: -set_initial_cache! = make_set_initial_cache(soil) -set_initial_cache!(p, Y, t0); - -# Timestepping functions: -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil) -jacobian! = ClimaLand.make_jacobian(soil) -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!) - timestepper = CTS.ARS111() ode_algo = CTS.IMEXAlgorithm( timestepper, @@ -176,31 +162,28 @@ ode_algo = CTS.IMEXAlgorithm( ), ) -# Define the problem and callbacks: -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -) -saveat = Array(t0:3600.0:tf) +saveat = Array(start_date:Hour(1):end_date); sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) updateat = deepcopy(saveat) -model_drivers = ClimaLand.get_drivers(soil) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); +simulation = LandSimulation( + start_date, + end_date, + dt, + soil; + set_ic! = set_ic!, + updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (saving_cb,), + diagnostics = (), +); # Solve -sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); +sol = solve!(simulation); # Figures @@ -223,14 +206,14 @@ ax = Axis( ) CairoMakie.lines!( ax, - sol.t ./ 3600 ./ 24, + FT.(sol.t) ./ 3600 ./ 24, sub .* (1000 * 3600 * 24), label = "Sublimation", color = :blue, ) CairoMakie.lines!( ax, - sol.t ./ 3600 ./ 24, + FT.(sol.t) ./ 3600 ./ 24, evap .* (1000 * 3600 * 24), label = "Evaporation", color = :black, diff --git a/docs/src/tutorials/standalone/Usage/model_tutorial.jl b/docs/src/tutorials/standalone/Usage/model_tutorial.jl index d160ffc652..4f9dafb725 100644 --- a/docs/src/tutorials/standalone/Usage/model_tutorial.jl +++ b/docs/src/tutorials/standalone/Usage/model_tutorial.jl @@ -124,6 +124,7 @@ using SciMLBase using Plots using ClimaCore using ClimaLand +import ClimaLand.Simulations: LandSimulation, solve! # Import the functions we are extending for our model: import ClimaLand: @@ -145,7 +146,7 @@ import ClimaLand: # the surface flux boundary value, the floating point precision, # and the domain of the model # (single column, global run, etc..). -struct RichardsTutorialModel{FT, D} <: AbstractModel{FT} +struct RichardsTutorialModel{FT, D} <: AbstractExpModel{FT} "van Genuchten model parameters" vGmodel::ClimaLand.Soil.vanGenuchten{FT} "Porosity [unitless]" @@ -322,34 +323,11 @@ soil = RichardsTutorialModel{Float32, typeof(domain)}( domain, ); -# Create the initial state structure, using the default method. This step -# creates the vector Y and cache p, but initializes them with zeros. -Y, p, cds = initialize(soil); - -# Note that `Y` has the structure we planned on in our `compute_exp_tendency!` function, for `x`, -Y.soil - -# The same is true for `p`: -p.soil - -# Here we now update `Y` in place with initial conditions of our choosing. - -Y.soil.θ = 0.25f0; - -# Set initial cache variable values, and inspect values: - -set_initial_cache! = make_set_initial_cache(soil); -set_initial_cache!(p, Y, 0.0); -@show p.soil.K - -@show p.soil.ψ - -@show p.soil.top_flux - -# Next up is to create the `exp_tendency!` function: -exp_tendency! = make_exp_tendency(soil); - -# # Running the simulation +# Here we create a function to update `Y` in place with initial conditions of our choosing. +function set_ic!(Y, p, t0, model) + Y.soil.θ = 0.25f0 +end +# # Creating and running the simulation # Set the initial and end times, timestep: @@ -361,16 +339,34 @@ dt = 1800.0; timestepper = CTS.RK4() ode_algo = CTS.ExplicitAlgorithm(timestepper) -# SciMLBase problem statement using CTS.jl internals: -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction(T_exp! = exp_tendency!), - Y, - (t0, tf), - p, +simulation = LandSimulation( + t0, + tf, + dt, + soil; + set_ic! = set_ic!, + solver_kwargs = (; save_everystep = true), + timestepper = ode_algo, + user_callbacks = (), + diagnostics = (), ); +Y = simulation._integrator.u +p = simulation._integrator.p +# Note that `Y` has the structure we planned on in our `compute_exp_tendency!` function, for `x`, +Y.soil + +# The same is true for `p`: +p.soil + +@show p.soil.K + +@show p.soil.ψ + +@show p.soil.top_flux + # Solve command: -sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = collect(t0:dt:tf)); +sol = solve!(simulation); # The solution is stored in `sol.u[k].soil.θ`, where k ranges over the number # of timesteps. diff --git a/experiments/calibration/model_interface.jl b/experiments/calibration/model_interface.jl index 0c091e0aec..251345cf88 100644 --- a/experiments/calibration/model_interface.jl +++ b/experiments/calibration/model_interface.jl @@ -359,7 +359,7 @@ function ClimaCalibrate.forward_model(iteration, member) calibrate_param_dict, ) - simulation = LandSimulation(FT, start_date, stop_date, Δt, model; outdir) + simulation = LandSimulation(start_date, stop_date, Δt, model; outdir) @info "Run: Global Soil-Canopy-Snow Model" @info "Resolution: $nelements" @info "Timestep: $Δt s" diff --git a/experiments/integrated/fluxnet/ozark_pft.jl b/experiments/integrated/fluxnet/ozark_pft.jl index c747d36273..28fd3488e5 100644 --- a/experiments/integrated/fluxnet/ozark_pft.jl +++ b/experiments/integrated/fluxnet/ozark_pft.jl @@ -19,6 +19,7 @@ using ClimaLand.Soil using ClimaLand.Soil.Biogeochemistry using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand import ClimaLand.Parameters as LP using ClimaDiagnostics @@ -105,12 +106,13 @@ canopy_domain = ClimaLand.Domains.obtain_surface_domain(land_domain) prognostic_land_components = (:canopy, :soil, :soilco2) # Set up the timestepping information for the simulation -t0 = Float64(0) dt = Float64(450) # 7.5 minutes N_spinup_days = 15 N_days = N_spinup_days + 340 -tf = Float64(t0 + N_days * 3600 * 24) -t_spinup = Float64(t0 + N_spinup_days * 3600 * 24) +N_seconds = Float64(N_days * 3600 * 24) +start_date = DateTime(2010) + Hour(time_offset) +end_date = start_date + Day(N_days) + timestepper = CTS.ARS111(); ode_algo = CTS.IMEXAlgorithm( @@ -164,7 +166,7 @@ pft_pcts = [ ac_canopy = ac_canopy * 3 # This reads in the data from the flux tower site and creates # the atmospheric and radiative driver structs for the model -start_date = DateTime(2010) + Hour(time_offset) + (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet( site_ID, lat, @@ -244,8 +246,8 @@ photosynthesis = FarquharModel{FT}(canopy_domain; photosynthesis_parameters) # Read in LAI from MODIS data surface_space = land_domain.space.surface; modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths( - start_date = start_date + Second(t0), - end_date = start_date + Second(t0) + Second(tf); + start_date = start_date, + end_date = end_date; context = ClimaComms.context(surface_space), ); LAI = ClimaLand.prescribed_lai_modis( @@ -297,16 +299,17 @@ canopy = Canopy.CanopyModel{FT}( # Integrated plant hydraulics and soil model land = SoilCanopyModel{FT}(soilco2, soil, canopy) -Y, p, cds = initialize(land) -exp_tendency! = make_exp_tendency(land) -imp_tendency! = make_imp_tendency(land); -jacobian! = make_jacobian(land); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); -FluxnetSimulations.set_fluxnet_ic!(Y, site_ID, start_date, time_offset, land) -set_initial_cache! = make_set_initial_cache(land) -set_initial_cache!(p, Y, t0); + + + +set_ic!(Y, _, _, model) = FluxnetSimulations.set_fluxnet_ic!( + Y, + site_ID, + start_date, + time_offset, + model, +) # Callbacks output_writer = ClimaDiagnostics.Writers.DictWriter() @@ -337,33 +340,20 @@ diags = ClimaLand.default_diagnostics( average_period = :hourly, ) -diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0, dt = dt); - -diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler); - -## How often we want to update the drivers. Note that this uses the defined `t0` and `tf` +## How often we want to update the drivers ## defined in the simulatons file data_dt = Float64(FluxnetSimulations.get_data_dt(site_ID)) -updateat = Array(t0:data_dt:tf) -model_drivers = ClimaLand.get_drivers(land) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, diag_cb) - - -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); - -sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb); +updateat = Array(start_date:Second(data_dt):end_date) +simulation = LandSimulation( + start_date, + end_date, + dt, + land; + set_ic! = set_ic!, + updateat, + diagnostics = diags, +) +sol = solve!(simulation) ClimaLand.Diagnostics.close_output_writers(diags) comparison_data = FluxnetSimulations.get_comparison_data(site_ID, time_offset) diff --git a/experiments/integrated/fluxnet/ozark_soilsnow.jl b/experiments/integrated/fluxnet/ozark_soilsnow.jl index 62cc6fea31..9b34c8ca94 100644 --- a/experiments/integrated/fluxnet/ozark_soilsnow.jl +++ b/experiments/integrated/fluxnet/ozark_soilsnow.jl @@ -18,6 +18,7 @@ using ClimaLand.Domains: Column using ClimaLand.Soil using ClimaLand.Snow import ClimaLand +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand.Parameters as LP import ClimaParams @@ -95,14 +96,13 @@ compartment_surfaces = n_stem > 0 ? [zmax, h_stem, h_canopy] : [zmax, h_leaf] domain = Column(; zlim = (zmin, zmax), nelements = nelements, dz_tuple) # Set up the timestepping information for the simulation -t0 = FT(1800) +start_date = DateTime(2010) + Hour(time_offset) N_days = 360 +end_date = start_date + Day(N_days) dt = FT(900) -tf = t0 + FT(3600 * 24 * N_days) -# This reads in the data from the flux tower site and creates -# the atmospheric and radiative driver structs for the model -start_date = DateTime(2010) + Hour(time_offset) +# Height of sensor on flux tower +atmos_h = FT(32) forcing = FluxnetSimulations.prescribed_forcing_fluxnet( site_ID, lat, @@ -159,61 +159,39 @@ snow_model = Snow.SnowModel( # Construct the land model land = ClimaLand.SoilSnowModel{FT}(; snow = snow_model, soil = soil_model) -Y, p, cds = initialize(land) -exp_tendency! = make_exp_tendency(land) -imp_tendency! = make_imp_tendency(land) -jacobian! = ClimaLand.make_jacobian(land); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); -set_initial_cache! = make_set_initial_cache(land) # Initial conditions -FluxnetSimulations.set_fluxnet_ic!(Y, site_ID, start_date, time_offset, land) -set_initial_cache!(p, Y, t0) +set_ic!(Y, p, _, model) = FluxnetSimulations.set_fluxnet_ic!( + Y, + site_ID, + start_date, + time_offset, + model, +) -saveat = Array(t0:dt:tf) +saveat = Array(start_date:Second(dt):end_date) sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) updateat = deepcopy(saveat) -model_drivers = ClimaLand.get_drivers(land) -updatefunc = ClimaLand.make_update_drivers(model_drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb) -# TIME STEPPING -timestepper = CTS.ARS111(); -ode_algo = CTS.IMEXAlgorithm( - timestepper, - CTS.NewtonsMethod( - max_iters = 3, - update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), - ), -); - -# Problem definition and callbacks -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); -sol = SciMLBase.solve( - prob, - ode_algo; - callback = cb, - dt = dt, - saveat = collect(t0:dt:tf), -); + +simulation = LandSimulation( + start_date, + end_date, + dt, + land; + user_callbacks = (saving_cb,), + set_ic! = set_ic!, + updateat, + solver_kwargs = (; saveat = saveat), +) +sol = solve!(simulation) # Plotting -daily = sol.t ./ 3600 ./ 24 +daily = FT.(sol.t) ./ 3600 ./ 24 savedir = joinpath(climaland_dir, "experiments/integrated/fluxnet/ozark_soilsnow") mkpath(savedir) @@ -227,7 +205,7 @@ ax1 = Axis(fig[2, 2], ylabel = "SWC", xlabel = "Days") lines!( ax1, daily, - [parent(sol.u[k].soil.ϑ_l)[end - 2] for k in 1:1:length(sol.t)], + [parent(sol.u[k].soil.ϑ_l)[end - 2] for k in 1:1:length(FT.(sol.t))], label = "10cm", ) lines!( @@ -235,7 +213,7 @@ lines!( daily, [ parent(sol.u[k].soil.θ_i .+ sol.u[k].soil.ϑ_l)[end - 2] for - k in 1:1:length(sol.t) + k in 1:1:length(FT.(sol.t)) ], label = "10cm, liq+ice", ) @@ -259,22 +237,23 @@ lines!( ) axislegend(ax2, position = :rb) ax3 = Axis(fig[2, 1], ylabel = "SWE (m)", xlabel = "Days") -lines!(ax3, daily, [parent(sol.u[k].snow.S)[1] for k in 1:1:length(sol.t)]) +lines!(ax3, daily, [parent(sol.u[k].snow.S)[1] for k in 1:1:length(FT.(sol.t))]) # Temp +sv_times = Dates.value.(Second.(sv.t .- sv.t[1])) ax4 = Axis(fig[1, 1], ylabel = "T (K)") hidexdecorations!(ax4, ticks = false) lines!( ax4, - sv.t ./ 24 ./ 3600, - [parent(sv.saveval[k].soil.T)[end - 2] for k in 1:1:length(sv.t)], + sv_times ./ 24 ./ 3600, + [parent(sv.saveval[k].soil.T)[end - 2] for k in 1:1:length(sv_times)], label = "Model 10 cm", ) lines!( ax4, - sv.t ./ 24 ./ 3600, - [parent(sv.saveval[k].snow.T)[1] for k in 1:1:length(sv.t)], + sv_times ./ 24 ./ 3600, + [parent(sv.saveval[k].snow.T)[1] for k in 1:1:length(sv_times)], label = "Snow", ) lines!( @@ -355,9 +334,9 @@ end compute_atmos_energy_fluxes(sv.saveval[k]) .- compute_energy_of_runoff(sv.saveval[k]) .- sv.saveval[k].soil.bottom_bc.heat, - )[end] for k in 1:1:(length(sv.t) - 1) + )[end] for k in 1:1:(length(sv_times) - 1) ], - ) * (sv.t[2] - sv.t[1]) + ) * (sv_times[2] - sv_times[1]) E_measured = [ sum(sol.u[k].soil.ρe_int) + parent(sol.u[k].snow.U)[end] for k in 1:1:length(sv.t) @@ -369,13 +348,13 @@ E_measured = [ compute_atmos_water_vol_fluxes(sv.saveval[k]) .- compute_runoff(sv.saveval[k]) .- sv.saveval[k].soil.bottom_bc.water, - )[end] for k in 1:1:(length(sv.t) - 1) + )[end] for k in 1:1:(length(sv_times) - 1) ], - ) * (sv.t[2] - sv.t[1]) + ) * (sv_times[2] - sv_times[1]) W_measured = [ sum(sol.u[k].soil.ϑ_l) + sum(sol.u[k].soil.θ_i) * _ρ_i / _ρ_l + - parent(sol.u[k].snow.S)[end] for k in 1:1:length(sv.t) + parent(sol.u[k].snow.S)[end] for k in 1:1:length(sv_times) ] lines!( ax1, diff --git a/experiments/integrated/fluxnet/run_fluxnet.jl b/experiments/integrated/fluxnet/run_fluxnet.jl index dc295e9db1..6960bca4d5 100644 --- a/experiments/integrated/fluxnet/run_fluxnet.jl +++ b/experiments/integrated/fluxnet/run_fluxnet.jl @@ -17,6 +17,7 @@ using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics import ClimaLand import ClimaLand.Parameters as LP +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand.FluxnetSimulations as FluxnetSimulations using CairoMakie, ClimaAnalysis, GeoMakie, Poppler_jll, Printf, StatsBase @@ -100,12 +101,9 @@ land_domain = Column(; surface_domain = ClimaLand.Domains.obtain_surface_domain(land_domain) # Set up the timestepping information for the simulation -t0 = Float64(0) dt = Float64(450) # 7.5 minutes N_spinup_days = 15 N_days = N_spinup_days + 340 -tf = Float64(t0 + N_days * 3600 * 24) -t_spinup = Float64(t0 + N_spinup_days * 3600 * 24) timestepper = CTS.ARS111(); ode_algo = CTS.IMEXAlgorithm( @@ -118,7 +116,11 @@ ode_algo = CTS.IMEXAlgorithm( # This reads in the data from the flux tower site and creates # the atmospheric and radiative driver structs for the model -(start_date, end_date) = FluxnetSimulations.get_data_dates(site_ID, time_offset) +(start_date, end_date) = FluxnetSimulations.get_data_dates( + site_ID, + time_offset; + duration = Day(N_days), +) (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet( site_ID, lat, @@ -200,8 +202,8 @@ photosynthesis = FarquharModel{FT}(surface_domain; photosynthesis_parameters) surface_space = land_domain.space.surface; modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(; context = ClimaComms.context(surface_space), - start_date = start_date + Second(t0), - end_date = start_date + Second(t0) + Second(tf), + start_date, + end_date, ) LAI = ClimaLand.prescribed_lai_modis( modis_lai_ncdata_path, @@ -263,17 +265,15 @@ snow = Snow.SnowModel( # Integrated plant hydraulics, soil, and snow model land = LandModel{FT}(canopy, snow, soil, soilco2); -Y, p, cds = initialize(land); -FluxnetSimulations.set_fluxnet_ic!(Y, site_ID, start_date, time_offset, land) -set_initial_cache! = make_set_initial_cache(land); -set_initial_cache!(p, Y, t0); +set_ic!(Y, p, _, model) = FluxnetSimulations.set_fluxnet_ic!( + Y, + site_ID, + start_date, + time_offset, + model, +) -exp_tendency! = make_exp_tendency(land); -imp_tendency! = make_imp_tendency(land); -jacobian! = make_jacobian(land); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); # Callbacks @@ -304,32 +304,20 @@ diags = ClimaLand.default_diagnostics( average_period = :hourly, ); -diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0, dt = dt); - -diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler); - -## How often we want to update the drivers. Note that this uses the defined `t0`, and `tf` -## defined in the simulatons file -data_dt = Float64(FluxnetSimulations.get_data_dt(site_ID)); -updateat = Array(t0:data_dt:tf); -model_drivers = ClimaLand.get_drivers(land); -updatefunc = ClimaLand.make_update_drivers(model_drivers); -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc); -cb = SciMLBase.CallbackSet(driver_cb, diag_cb); - -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); +## How often we want to update the drivers. +data_dt = Second(FluxnetSimulations.get_data_dt(site_ID)) -@time sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb); +updateat = Array(start_date:data_dt:end_date) +simulation = LandSimulation( + start_date, + end_date, + dt, + land; + set_ic! = set_ic!, + updateat, + diagnostics = diags, +) +@time sol = solve!(simulation) ClimaLand.Diagnostics.close_output_writers(diags) comparison_data = FluxnetSimulations.get_comparison_data(site_ID, time_offset) diff --git a/experiments/integrated/global/global_soil_canopy.jl b/experiments/integrated/global/global_soil_canopy.jl index 5c0bc0c0c1..c165c2e8f8 100644 --- a/experiments/integrated/global/global_soil_canopy.jl +++ b/experiments/integrated/global/global_soil_canopy.jl @@ -15,6 +15,7 @@ using ClimaLand.Soil using ClimaLand.Canopy import ClimaLand import ClimaLand.Parameters as LP +import ClimaLand.Simulations: LandSimulation, solve! import CairoMakie import GeoMakie @@ -47,9 +48,8 @@ subsurface_space = domain.space.subsurface # Set up dates and times for the simulation start_date = DateTime(2008); -t0 = 0.0 dt = 450.0 -tf = 3600 +end_date = start_date + Dates.Second(3600) # Forcing data era5_ncdata_path = @@ -91,8 +91,8 @@ canopy_forcing = (; atmos, radiation, ground) # Set up plant hydraulics modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(; context = nothing, - start_date = start_date + Second(t0), - end_date = start_date + Second(t0) + Second(tf), + start_date = start_date, + end_date = end_date, ) LAI = ClimaLand.prescribed_lai_modis( modis_lai_ncdata_path, @@ -111,19 +111,10 @@ canopy = Canopy.CanopyModel{FT}( # Combine the soil and canopy models into a single prognostic land model land = SoilCanopyModel{FT}(soilco2, soil, canopy) -Y, p, cds = initialize(land) - ic_path = ClimaLand.Artifacts.soil_ic_2008_50m_path(; context = context) -set_initial_state! = - ClimaLand.Simulations.make_set_initial_state_from_file(ic_path, land) -set_initial_cache! = make_set_initial_cache(land) -exp_tendency! = make_exp_tendency(land) -imp_tendency! = ClimaLand.make_imp_tendency(land) -jacobian! = ClimaLand.make_jacobian(land) - -set_initial_state!(Y, p, t0, land) -set_initial_cache!(p, Y, t0) - +set_ic! = + (Y, p, t0, model) -> + ClimaLand.Simulations.make_set_initial_state_from_file(ic_path, model) stepper = CTS.ARS343() ode_algo = CTS.IMEXAlgorithm( stepper, @@ -133,21 +124,6 @@ ode_algo = CTS.IMEXAlgorithm( ), ) -# set up jacobian info -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!) - -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -) - # ClimaDiagnostics nc_writer = ClimaDiagnostics.Writers.NetCDFWriter(subsurface_space, outdir; start_date) @@ -159,23 +135,17 @@ diags = ClimaLand.default_diagnostics( average_period = :hourly, ) -diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = dt) - -diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) - -updateat = Array(t0:(3600 * 3):tf) -drivers = ClimaLand.get_drivers(land) -updatefunc = ClimaLand.make_update_drivers(drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) - -sol = ClimaComms.@time ClimaComms.device() SciMLBase.solve( - prob, - ode_algo; - dt = dt, - callback = SciMLBase.CallbackSet(driver_cb, diag_cb), +simulation = LandSimulation( + start_date, + end_date, + dt, + land; + outdir, + diagnostics = diags, + timestepper = ode_algo, + user_callbacks = (), ) - +ClimaLand.Simulations.solve!(simulation) ClimaLand.Diagnostics.close_output_writers(diags) # ClimaAnalysis diff --git a/experiments/integrated/performance/integrated_timestep_test.jl b/experiments/integrated/performance/integrated_timestep_test.jl index 4a46b1df1f..1227d6ff58 100644 --- a/experiments/integrated/performance/integrated_timestep_test.jl +++ b/experiments/integrated/performance/integrated_timestep_test.jl @@ -47,7 +47,7 @@ import ClimaComms ClimaComms.@import_required_backends import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput import ClimaUtilities.OutputPathGenerator: generate_output_path - +import ClimaUtilities.TimeManager: date using ClimaLand using ClimaLand.Domains using ClimaLand.Soil @@ -56,6 +56,7 @@ using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics import ClimaLand import ClimaLand.Parameters as LP +import ClimaLand.Simulations: LandSimulation, solve! import ClimaParams @@ -67,8 +68,7 @@ try catch end -function set_initial_conditions(land, t0) - Y, p, cds = initialize(land) +function set_ic!(Y, p, t0, model) FT = eltype(Y.soil.ϑ_l) set_initial_cache! = make_set_initial_cache(land) @@ -110,8 +110,7 @@ function set_initial_conditions(land, t0) end Y.canopy.energy.T = FT(297.5) - set_initial_cache!(p, Y, t0) - return Y, p + return end context = ClimaComms.context() @@ -303,14 +302,10 @@ canopy = Canopy.CanopyModel{FT}( # Integrated plant and soil model land = SoilCanopyModel{FT}(soilco2, soil, canopy) -# Define explicit and implicit tendencies, and the jacobian -exp_tendency! = make_exp_tendency(land) -imp_tendency! = make_imp_tendency(land); -jacobian! = make_jacobian(land); - # Timestepping information N_hours = 8 tf = Float64(t0 + N_hours * 3600.0) +stop_date = start_date + Dates.Second(round(tf)) sim_time = round((tf - t0) / 3600, digits = 2) # simulation length in hours timestepper = CTS.ARS111() @@ -322,10 +317,6 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Set up simulation callbacks -drivers = ClimaLand.get_drivers(land) -updatefunc = ClimaLand.make_update_drivers(drivers) - # Choose timesteps and set up arrays to store information ref_dt = 50.0 dts = [ref_dt, 100.0, 200.0, 450.0, 900.0, 1800.0, 3600.0] @@ -340,36 +331,23 @@ times = [] for dt in dts @info dt # Initialize model and set initial conditions before each simulation - Y, p = set_initial_conditions(land, t0) - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, - ) - - # Create callback for driver updates - saveat = vcat(Array(t0:(3 * 3600):tf), tf) - updateat = vcat(Array(t0:(3 * 3600):tf), tf) - driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) - - # Solve simulation - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, + simulation = LandSimulation( + start_date, + stop_date, + dt, + land; + solver_kwargs = (; + saveat = vcat( + collect(start_date:Dates.Hour(3):stop_date), + stop_date, + ) ), - Y, - (t0, tf), - p, - ) - @time sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - callback = driver_cb, - saveat = saveat, + diagnostics = [], + updateat = vcat(collect(start_date:Dates.Hour(3):stop_date), stop_date), + timestepper = ode_algo, + set_ic!, ) - + sol = ClimaLand.Simulations.solve!(simulation) # Save results for comparison if dt == ref_dt global ref_T = @@ -430,9 +408,8 @@ ax3 = Axis( ylabel = "T (K)", title = "T throughout simulation; length = $(sim_time) hours, 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") + lines!(ax3, FT.(times[i]) ./ 3600.0, T_states[i], label = "dt $(dts[i]) s") end axislegend(ax3, position = :rt) save(joinpath(savedir, "states.png"), fig3) diff --git a/experiments/long_runs/bucket.jl b/experiments/long_runs/bucket.jl index 33a9c2c94e..003239bac3 100644 --- a/experiments/long_runs/bucket.jl +++ b/experiments/long_runs/bucket.jl @@ -138,7 +138,6 @@ timestepper = CTS.ExplicitAlgorithm(timestepper) # Create the simulation simulation = LandSimulation( - FT, start_date, stop_date, Δt, diff --git a/experiments/long_runs/land_region.jl b/experiments/long_runs/land_region.jl index 4614b06749..10c754f8c3 100644 --- a/experiments/long_runs/land_region.jl +++ b/experiments/long_runs/land_region.jl @@ -163,7 +163,7 @@ domain = ClimaLand.Domains.HybridBox(; ) params = LP.LandParameters(FT) model = setup_model(FT, context, start_date, Δt, domain, params) -simulation = LandSimulation(FT, start_date, stop_date, Δt, model; outdir) +simulation = LandSimulation(start_date, stop_date, Δt, model; outdir) ClimaLand.Simulations.solve!(simulation) LandSimVis.make_heatmaps(simulation; savedir = root_path, date = stop_date) diff --git a/experiments/long_runs/low_res_snowy_land.jl b/experiments/long_runs/low_res_snowy_land.jl index 6043e57d0e..b07bdaab9e 100644 --- a/experiments/long_runs/low_res_snowy_land.jl +++ b/experiments/long_runs/low_res_snowy_land.jl @@ -172,7 +172,7 @@ user_callbacks = ( ClimaLand.ReportCallback(10000), ) simulation = - LandSimulation(FT, start_date, stop_date, Δt, model; user_callbacks, outdir) + LandSimulation(start_date, stop_date, Δt, model; user_callbacks, outdir) @info "Run: Global Soil-Canopy-Snow Model" @info "Resolution: $nelements" @info "Timestep: $Δt s" diff --git a/experiments/long_runs/snowy_land.jl b/experiments/long_runs/snowy_land.jl index 264646614d..d7cea75816 100644 --- a/experiments/long_runs/snowy_land.jl +++ b/experiments/long_runs/snowy_land.jl @@ -189,7 +189,7 @@ domain = ClimaLand.Domains.global_domain( ) params = LP.LandParameters(FT) model = setup_model(FT, start_date, stop_date, Δt, domain, params) -simulation = LandSimulation(FT, start_date, stop_date, Δt, model; outdir) +simulation = LandSimulation(start_date, stop_date, Δt, model; outdir) @info "Run: Global Soil-Canopy-Snow Model" @info "Resolution: $nelements" @info "Timestep: $Δt s" diff --git a/experiments/long_runs/soil.jl b/experiments/long_runs/soil.jl index eb2c0a374a..50f37e1039 100644 --- a/experiments/long_runs/soil.jl +++ b/experiments/long_runs/soil.jl @@ -93,7 +93,7 @@ diagnostics = ClimaLand.default_diagnostics( conservation_period = Day(10), ) simulation = - LandSimulation(FT, start_date, stop_date, Δt, model; outdir, diagnostics) + LandSimulation(start_date, stop_date, Δt, model; outdir, diagnostics) @info "Run: Global Soil Model" @info "Resolution: $nelements" diff --git a/experiments/standalone/Biogeochemistry/experiment.jl b/experiments/standalone/Biogeochemistry/experiment.jl index 8d36a9f43e..a11af2f1c1 100644 --- a/experiments/standalone/Biogeochemistry/experiment.jl +++ b/experiments/standalone/Biogeochemistry/experiment.jl @@ -8,6 +8,7 @@ using ClimaLand.Domains: Column using ClimaLand.Soil using ClimaLand.Soil.Biogeochemistry using ClimaLand.Soil.Biogeochemistry: MicrobeProduction +import ClimaLand.Simulations: LandSimulation, solve! using Dates import ClimaLand.Parameters as LP @@ -117,9 +118,6 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf)) soilco2_args = soilco2_args, ) - Y, p, coords = initialize(model) - set_initial_cache! = make_set_initial_cache(model) - function init_soil!(Y, z, params) ν = params.ν FT = eltype(Y.soil.ϑ_l) @@ -149,11 +147,11 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf)) Y.soilco2.C .= CO2_profile.(z) end - z = coords.subsurface.z - init_soil!(Y, z, model.soil.parameters) - init_co2!(Y, z) - set_initial_cache!(p, Y, t0) - Soil_bio_exp_tendency! = make_exp_tendency(model) + function set_ic!(Y, p, t0, model) + z = model.soil.domain.fields.z + init_co2!(Y, z) + init_soil!(Y, z, model.soil.parameters) + end timestepper = CTS.RK4() ode_algo = CTS.ExplicitAlgorithm(timestepper) @@ -165,18 +163,20 @@ for (FT, tf) in ((Float32, 2 * dt), (Float64, tf)) ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) updateat = deepcopy(saveat) - drivers = ClimaLand.get_drivers(model) - updatefunc = ClimaLand.make_update_drivers(drivers) - driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) - cb = SciMLBase.CallbackSet(driver_cb, saving_cb) - - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction(T_exp! = Soil_bio_exp_tendency!), - Y, - (t0, tf), - p, + + simulation = LandSimulation( + t0, + tf, + dt, + model; + diagnostics = [], + timestepper = ode_algo, + set_ic!, + user_callbacks = (saving_cb,), + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), ) - sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb) + sol = ClimaLand.Simulations.solve!(simulation) # Check that simulation still has correct float type @assert eltype(sol.u[end].soil) == FT diff --git a/experiments/standalone/Bucket/bucket_era5.jl b/experiments/standalone/Bucket/bucket_era5.jl index c1ad94b894..7d17578622 100644 --- a/experiments/standalone/Bucket/bucket_era5.jl +++ b/experiments/standalone/Bucket/bucket_era5.jl @@ -44,6 +44,7 @@ import ClimaLand.Parameters as LP using ClimaLand.Bucket: BucketModel, BucketModelParameters, PrescribedBaregroundAlbedo using ClimaLand.Domains: coordinates, Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand: initialize, make_update_aux, @@ -156,37 +157,22 @@ model = BucketModel( radiation = bucket_rad, ); -Y, p, coords = initialize(model); - -Y.bucket.T .= FT(270); -Y.bucket.W .= FT(0.05); -Y.bucket.Ws .= FT(0.0); -Y.bucket.σS .= FT(0.08); - -set_initial_cache! = make_set_initial_cache(model); -set_initial_cache!(p, Y, t0); -exp_tendency! = make_exp_tendency(model); -timestepper = CTS.RK4() -ode_algo = CTS.ExplicitAlgorithm(timestepper) -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!), - Y, - (t0, tf), - p, -); - -saveat = [promote(t0:(3 * Δt):tf...)...]; +saveat = [promote(t0:(3 * Δt):tf...)...] saved_values = (; t = Array{typeof(t0)}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ); saving_cb = ClimaLand.NonInterpSavingCallback(saved_values, saveat); -updateat = copy(saveat) -drivers = ClimaLand.get_drivers(model) -updatefunc = ClimaLand.make_update_drivers(drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -# Diagnostics +function set_ic!(Y, p, t0, model) + Y.bucket.T .= FT(270) + Y.bucket.W .= FT(0.05) + Y.bucket.Ws .= FT(0.0) + Y.bucket.σS .= FT(0.08) +end + +timestepper = CTS.RK4() +ode_algo = CTS.ExplicitAlgorithm(timestepper) nc_writer = ClimaDiagnostics.Writers.NetCDFWriter( subsurface_space, output_dir; @@ -197,20 +183,25 @@ diags = ClimaLand.default_diagnostics( start_date; output_writer = nc_writer, average_period = :daily, -) +); -diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) +simulation = LandSimulation( + t0, + tf, + Δt, + model; + outdir = output_dir, + solver_kwargs = (; saveat = copy(saveat)), + updateat = copy(saveat), + user_callbacks = (saving_cb,), + set_ic!, + timestepper = ode_algo, + diagnostics = diags, +) -diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb, diag_cb) -sol = ClimaComms.@time ClimaComms.device() SciMLBase.solve( - prob, - ode_algo; - dt = Δt, - saveat = saveat, - callback = cb, +sol = ClimaComms.@time ClimaComms.device() ClimaLand.Simulations.solve!( + simulation, ); ClimaLand.Diagnostics.close_output_writers(diags) @@ -233,7 +224,7 @@ for short_name in short_names end # Interpolate to grid -space = axes(coords.surface) +space = simulation.model.domain.space.surface num_pts = 21 if regional_simulation radius_earth = FT(6.378e6) @@ -272,7 +263,9 @@ remapper = Remapping.Remapper(space, hcoords) W = Array(Remapping.interpolate(remapper, sol.u[end].bucket.W)) Ws = Array(Remapping.interpolate(remapper, sol.u[end].bucket.Ws)) σS = Array(Remapping.interpolate(remapper, sol.u[end].bucket.σS)) -T_sfc = Array(Remapping.interpolate(remapper, prob.p.bucket.T_sfc)) +T_sfc = Array( + Remapping.interpolate(remapper, simulation._integrator.p.bucket.T_sfc), +) # save prognostic state to CSV - for comparison between GPU and CPU output open(joinpath(output_dir, "tf_state_$(device_suffix)_era5.txt"), "w") do io diff --git a/experiments/standalone/Bucket/global_bucket_function.jl b/experiments/standalone/Bucket/global_bucket_function.jl index 2580dc5c36..63ec489c31 100644 --- a/experiments/standalone/Bucket/global_bucket_function.jl +++ b/experiments/standalone/Bucket/global_bucket_function.jl @@ -34,6 +34,7 @@ import ClimaParams using ClimaLand.Bucket: BucketModel, BucketModelParameters, PrescribedBaregroundAlbedo using ClimaLand.Domains: coordinates, Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand: initialize, make_update_aux, @@ -92,18 +93,17 @@ z_0b = FT(1e-3); κ_soil = FT(0.7); ρc_soil = FT(2e6); τc = FT(3600); -t0 = 0.0; -tf = 7 * 86400; Δt = 3600.0; bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc); start_date = DateTime(2005); +end_date = start_date + Second(7 * 86400); # Precipitation: precip = (t) -> 0; -snow_precip = (t) -> -5e-7 * (t < 1 * 86400); +snow_precip = (t) -> -5e-7 * (FT(t) < 1 * 86400); # Diurnal temperature variations: -T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2); +T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * FT(t) / 86400 - π / 2); # Constant otherwise: u_atmos = (t) -> 3.0; q_atmos = (t) -> 0.001; @@ -126,8 +126,8 @@ bucket_atmos = PrescribedAtmosphere( # peak at local noon, and a prescribed downwelling LW radiative # flux, assuming the air temperature is on average 275 degrees # K with a diurnal amplitude of 5 degrees K: -SW_d = (t) -> max(1361 * sin(2π * t / 86400 - π / 2), 0.0); -LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2))^4; +SW_d = (t) -> max(1361 * sin(2π * FT(t) / 86400 - π / 2), 0.0); +LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * FT(t) / 86400 - π / 2))^4; bucket_rad = PrescribedRadiativeFluxes( FT, TimeVaryingInput(SW_d), @@ -142,26 +142,6 @@ model = BucketModel( radiation = bucket_rad, ); -Y, p, coords = initialize(model); - -Y.bucket.T .= FT(270); -Y.bucket.W .= FT(0.05); -Y.bucket.Ws .= FT(0.0); -Y.bucket.σS .= FT(0.08); - -set_initial_cache! = make_set_initial_cache(model); -set_initial_cache!(p, Y, t0); -exp_tendency! = make_exp_tendency(model); -timestepper = CTS.RK4() -ode_algo = CTS.ExplicitAlgorithm(timestepper) -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction(T_exp! = exp_tendency!, dss! = ClimaLand.dss!), - Y, - (t0, tf), - p, -); - -# ClimaDiagnostics output_dir = ClimaUtilities.OutputPathGenerator.generate_output_path(outdir) space = bucket_domain.space.subsurface @@ -175,23 +155,35 @@ diags = ClimaLand.default_diagnostics( average_period = :daily, ) +timestepper = CTS.RK4() +ode_algo = CTS.ExplicitAlgorithm(timestepper) + +function set_ic!(Y, p, t0, model) + Y.bucket.T .= FT(270) + Y.bucket.W .= FT(0.05) + Y.bucket.Ws .= FT(0.0) + Y.bucket.σS .= FT(0.08) +end + +simulation = LandSimulation( + start_date, + end_date, + Δt, + model; + outdir = output_dir, + updateat = collect(start_date:Second(Δt * 3):end_date), + set_ic!, + timestepper = ode_algo, + diagnostics = diags, +) + -diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) -diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) -updateat = collect(t0:(Δt * 3):tf); -drivers = ClimaLand.get_drivers(model) -updatefunc = ClimaLand.make_update_drivers(drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -sol = ClimaComms.@time ClimaComms.device() SciMLBase.solve( - prob, - ode_algo; - dt = Δt, - callback = SciMLBase.CallbackSet(driver_cb, diag_cb), -) +sol = ClimaComms.@time ClimaComms.device() ClimaLand.Simulations.solve!( + simulation, +); ClimaLand.Diagnostics.close_output_writers(diags) @@ -223,7 +215,7 @@ for short_name in vcat(short_names_2D..., short_names_3D...) end # Interpolate to grid -space = axes(coords.surface) +space = axes(sol.prob.p.drivers.T) longpts = range(-180.0, 180.0, 21) latpts = range(-90.0, 90.0, 21) hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts] @@ -232,7 +224,7 @@ remapper = Remapping.Remapper(space, hcoords) W = Array(Remapping.interpolate(remapper, sol.u[end].bucket.W)) Ws = Array(Remapping.interpolate(remapper, sol.u[end].bucket.Ws)) σS = Array(Remapping.interpolate(remapper, sol.u[end].bucket.σS)) -T_sfc = Array(Remapping.interpolate(remapper, prob.p.bucket.T_sfc)) +T_sfc = Array(Remapping.interpolate(remapper, sol.prob.p.bucket.T_sfc)) # save prognostic state to CSV - for comparison between diff --git a/experiments/standalone/Bucket/global_bucket_temporalmap.jl b/experiments/standalone/Bucket/global_bucket_temporalmap.jl index 8708c322e2..1ec1d3b594 100644 --- a/experiments/standalone/Bucket/global_bucket_temporalmap.jl +++ b/experiments/standalone/Bucket/global_bucket_temporalmap.jl @@ -40,6 +40,7 @@ import ClimaLand.Parameters as LP using ClimaLand.Bucket: BucketModel, BucketModelParameters, PrescribedSurfaceAlbedo using ClimaLand.Domains: coordinates, Column +import ClimaLand.Simulations: LandSimulation, solve! using ClimaLand: initialize, make_update_aux, @@ -78,12 +79,13 @@ device_suffix = typeof(context.device) <: ClimaComms.CPUSingleThreaded ? "cpu" : "gpu" outdir = "experiments/standalone/Bucket/artifacts_temporalmap_$(device_suffix)" t0 = 0.0; +start_date = DateTime(2005) +end_date = start_date + Second(50 * 86400) # run for 50 days to test monthly file update -tf = 50 * 86400; Δt = 3600.0; -function setup_prob(t0, tf, Δt, outdir) +function setup_prob(start_date, end_date, Δt, outdir) # We set up the problem in a function so that we can make multiple copies (for profiling) output_dir = ClimaUtilities.OutputPathGenerator.generate_output_path(outdir) @@ -96,7 +98,6 @@ function setup_prob(t0, tf, Δt, outdir) nelements = (10, 10), # this failed with (50,10) dz_tuple = FT.((1.0, 0.05)), ) - start_date = DateTime(2005) # Initialize parameters σS_c = FT(0.2) @@ -116,9 +117,9 @@ function setup_prob(t0, tf, Δt, outdir) # Precipitation: precip = (t) -> 0 - snow_precip = (t) -> -5e-7 * (t < 1 * 86400) + snow_precip = (t) -> -5e-7 * (FT(t) < 1 * 86400) # Diurnal temperature variations: - T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2) + T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * FT(t) / 86400 - π / 2) # Constant otherwise: u_atmos = (t) -> 3.0 q_atmos = (t) -> 0.001 @@ -141,8 +142,9 @@ function setup_prob(t0, tf, Δt, outdir) # peak at local noon, and a prescribed downwelling LW radiative # flux, assuming the air temperature is on average 275 degrees # K with a diurnal amplitude of 5 degrees K: - SW_d = (t) -> max(1361 * sin(2π * t / 86400 - π / 2), 0.0) - LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2))^4 + SW_d = (t) -> max(1361 * sin(2π * FT(t) / 86400 - π / 2), 0.0) + LW_d = + (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * FT(t) / 86400 - π / 2))^4 bucket_rad = PrescribedRadiativeFluxes( FT, TimeVaryingInput(SW_d), @@ -158,31 +160,23 @@ function setup_prob(t0, tf, Δt, outdir) radiation = bucket_rad, ) - Y, p, _coords = initialize(model) - - Y.bucket.T .= FT(270) - Y.bucket.W .= FT(0.05) - Y.bucket.Ws .= FT(0.0) - Y.bucket.σS .= FT(0.08) - - set_initial_cache! = make_set_initial_cache(model) - set_initial_cache!(p, Y, t0) - exp_tendency! = make_exp_tendency(model) - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction((T_exp!) = exp_tendency!, (dss!) = ClimaLand.dss!), - Y, - (t0, tf), - p, - ) - saveat = collect(t0:(Δt * 3):tf) + + function set_ic!(Y, p, t0, model) + Y.bucket.T .= FT(270) + Y.bucket.W .= FT(0.05) + Y.bucket.Ws .= FT(0.0) + Y.bucket.σS .= FT(0.08) + end + timestepper = CTS.RK4() + ode_algo = CTS.ExplicitAlgorithm(timestepper) + + + saveat = collect(start_date:Second(Δt * 3):end_date) saved_values = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(saved_values, saveat) - updateat = copy(saveat) - drivers = ClimaLand.get_drivers(model) - updatefunc = ClimaLand.make_update_drivers(drivers) nc_writer = ClimaDiagnostics.Writers.NetCDFWriter( subsurface_space, output_dir; @@ -195,27 +189,30 @@ function setup_prob(t0, tf, Δt, outdir) average_period = :daily, ) - diagnostic_handler = - ClimaDiagnostics.DiagnosticsHandler(diags, Y, p, t0; dt = Δt) - - diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) - driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) - cb = SciMLBase.CallbackSet(driver_cb, saving_cb, diag_cb) + simulation = LandSimulation( + start_date, + end_date, + Δt, + model; + outdir = output_dir, + updateat = copy(saveat), + solver_kwargs = (; saveat = deepcopy(saveat)), + user_callbacks = (saving_cb,), + set_ic!, + timestepper = ode_algo, + diagnostics = diags, + ) - return prob, cb, saveat, saved_values, nc_writer + return simulation, saved_values, nc_writer end -prob, cb, saveat, saved_values, nc_writer = setup_prob(t0, tf, Δt, outdir); -timestepper = CTS.RK4() -ode_algo = CTS.ExplicitAlgorithm(timestepper) - -sol = ClimaComms.@time ClimaComms.device() SciMLBase.solve( - prob, - ode_algo; - dt = Δt, - saveat = saveat, - callback = cb, -) +simulation, saved_values, nc_writer = + setup_prob(start_date, end_date, Δt, outdir); + + +sol = ClimaComms.@time ClimaComms.device() ClimaLand.Simulations.solve!( + simulation, +); close(nc_writer) output_dir = nc_writer.output_dir @@ -241,28 +238,20 @@ if PROFILING # We run only for one day for profiling tf = 86400.0 - prob, cb, saveat, _, nc_writer = setup_prob(t0, tf, Δt, outdir) - - Profile.@profile SciMLBase.solve( - prob, - ode_algo; - dt = Δt, - saveat = saveat, - callback = cb, - ) + simulation, saved_values, nc_writer = + setup_prob(start_date, start_date + Second(tf), Δt, outdir) + + Profile.@profile ClimaLand.Simulations.solve!(simulation) results = Profile.fetch() flame_file = joinpath(outdir, "flame_$device_suffix.html") ProfileCanvas.html_file(flame_file, results) @info "Save compute flame to $flame_file" close(nc_writer) - prob, cb, saveat, _, nc_writer = setup_prob(t0, tf, Δt, outdir) - Profile.Allocs.@profile sample_rate = 0.1 SciMLBase.solve( - prob, - ode_algo; - dt = Δt, - saveat = saveat, - callback = cb, + simulation, saved_values, nc_writer = + setup_prob(start_date, start_date + Second(tf), Δt, outdir) + Profile.Allocs.@profile sample_rate = 0.1 ClimaLand.Simulations.solve!( + simulation, ) results = Profile.Allocs.fetch() profile = ProfileCanvas.view_allocs(results) @@ -282,7 +271,7 @@ remapper = Remapping.Remapper(space, hcoords) W = Array(Remapping.interpolate(remapper, sol.u[end].bucket.W)) Ws = Array(Remapping.interpolate(remapper, sol.u[end].bucket.Ws)) σS = Array(Remapping.interpolate(remapper, sol.u[end].bucket.σS)) -T_sfc = Array(Remapping.interpolate(remapper, prob.p.bucket.T_sfc)) +T_sfc = Array(Remapping.interpolate(remapper, sol.prob.p.bucket.T_sfc)) # save prognostic state to CSV - for comparison between GPU and CPU output open( diff --git a/experiments/standalone/Soil/energy_hydrology_saturated.jl b/experiments/standalone/Soil/energy_hydrology_saturated.jl index 64b1eab9d6..6a2eb7ec7a 100644 --- a/experiments/standalone/Soil/energy_hydrology_saturated.jl +++ b/experiments/standalone/Soil/energy_hydrology_saturated.jl @@ -7,6 +7,7 @@ using ClimaLand using ClimaLand.Domains: Column using ClimaLand.Soil import ClimaLand +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand.Parameters as LP FT = Float32; @@ -61,14 +62,7 @@ soil = Soil.EnergyHydrology{FT}(; boundary_conditions = boundary_fluxes, sources = (), ) -Y, p, cds = initialize(soil); -set_initial_cache! = make_set_initial_cache(soil) -exp_tendency! = make_exp_tendency(soil) -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); -jac_kwargs = - (; jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), Wfact = jacobian!); function hydrostatic_equilibrium(z, z_interface, params) (; ν, S_s, hydrology_cm) = params @@ -79,7 +73,9 @@ function hydrostatic_equilibrium(z, z_interface, params) return ν * (1 + (α * (z - z_interface))^n)^(-m) end end -function init_soil!(Y, z, params) +function init_soil!(Y, p, t0, model) + params = model.parameters + z = model.domain.fields.z FT = eltype(Y.soil.ϑ_l) Y.soil.ϑ_l .= hydrostatic_equilibrium.(z, FT(-0.45), params) Y.soil.θ_i .= 0 @@ -106,29 +102,30 @@ ode_algo = CTS.IMEXAlgorithm( ), ); -# Define the problem and callbacks: -prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, -); +simulation_dt_small = + LandSimulation(t0, tf, FT(100), soil; set_ic! = init_soil!) + +simulation_dt_large = LandSimulation( + t0, + tf, + FT(900), + soil; + set_ic! = init_soil!, + timestepper = ode_algo, +) + # Solve at small dt -init_soil!(Y, z, soil.parameters); -set_initial_cache!(p, Y, t0); -sol_dt_small = SciMLBase.solve(prob, ode_algo; dt = 100.0); -init_soil!(Y, z, soil.parameters); -set_initial_cache!(p, Y, t0); -sol_dt_large = SciMLBase.solve(prob, ode_algo; dt = 900.0); -@assert sum(isnan.(sol_dt_large.u[end])) == 0 +sol_dt_small = solve!(simulation_dt_small); +sol_dt_large = solve!(simulation_dt_large); +@assert sum(isnan.(simulation_dt_large._integrator.u)) == 0 norm(x) = sqrt(mean(parent(x .^ 2))) rmse(x, y) = sqrt(mean(parent((x .- y) .^ 2))) -@assert rmse(sol_dt_small[end].soil.ϑ_l, sol_dt_large[end].soil.ϑ_l) / - norm(sol_dt_small[end].soil.ϑ_l) < sqrt(eps(FT)) -@assert rmse(sol_dt_small[end].soil.ρe_int, sol_dt_large[end].soil.ρe_int) / - norm(sol_dt_small[end].soil.ρe_int) < sqrt(eps(FT)) +@assert rmse( + simulation_dt_small._integrator.u.soil.ϑ_l, + simulation_dt_large._integrator.u.soil.ϑ_l, +) / norm(simulation_dt_small._integrator.u.soil.ϑ_l) < sqrt(eps(FT)) +@assert rmse( + simulation_dt_small._integrator.u.soil.ρe_int, + simulation_dt_large._integrator.u.soil.ρe_int, +) / norm(simulation_dt_small._integrator.u.soil.ρe_int) < sqrt(eps(FT)) diff --git a/experiments/standalone/Soil/richards_comparison.jl b/experiments/standalone/Soil/richards_comparison.jl index 3881ed7987..fcdb56fed0 100644 --- a/experiments/standalone/Soil/richards_comparison.jl +++ b/experiments/standalone/Soil/richards_comparison.jl @@ -12,6 +12,7 @@ using ClimaLand using ClimaLand.Domains: Column using ClimaLand.Soil import ClimaUtilities.OutputPathGenerator: generate_output_path +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand import ClimaLand.Parameters as LP @@ -67,16 +68,11 @@ outdir = generate_output_path( boundary_conditions = boundary_states, sources = sources, ) - set_initial_cache! = make_set_initial_cache(soil) - - Y, p, coords = initialize(soil) # specify ICs - Y.soil.ϑ_l .= FT(0.24) - exp_tendency! = make_exp_tendency(soil) - imp_tendency! = ClimaLand.make_imp_tendency(soil) - jacobian! = ClimaLand.make_jacobian(soil) - set_initial_cache!(p, Y, t0) + function set_ic!(Y, p, t0, model) + Y.soil.ϑ_l .= FT(0.24) + end stepper = CTS.ARS111() norm_condition = CTS.MaximumError(FT(1e-8)) @@ -90,28 +86,18 @@ outdir = generate_output_path( ), ) - # set up jacobian info - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, + simulation = LandSimulation( + t0, + tf, + dt, + soil; + outdir, + set_ic!, + timestepper = ode_algo, + solver_kwargs = (; saveat = collect(t0:10000:tf)), ) - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, - ) - sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - saveat = collect(t0:10000:tf), - ) + sol = solve!(simulation) # Check that simulation still has correct float type @assert eltype(sol.u[end].soil) == FT @@ -176,15 +162,11 @@ end sources = sources, ) - Y, p, coords = initialize(soil) - set_initial_cache! = make_set_initial_cache(soil) # specify ICs - Y.soil.ϑ_l .= FT(0.1) - exp_tendency! = make_exp_tendency(soil) - imp_tendency! = ClimaLand.make_imp_tendency(soil) - jacobian! = ClimaLand.make_jacobian(soil) - set_initial_cache!(p, Y, t0) + function set_ic!(Y, p, t0, model) + Y.soil.ϑ_l .= FT(0.1) + end stepper = CTS.ARS111() norm_condition = CTS.MaximumError(FT(1e-8)) @@ -197,28 +179,18 @@ end convergence_checker = conv_checker, ), ) - # set up jacobian info - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, + simulation = LandSimulation( + t0, + tf, + dt, + soil; + outdir, + set_ic!, + timestepper = ode_algo, + solver_kwargs = (; saveat = collect(t0:(60 * dt):tf)), ) - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, - ) - sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - saveat = collect(t0:(60 * dt):tf), - ) + sol = solve!(simulation) # Check that simulation still has correct float type @assert eltype(sol.u[end].soil) == FT diff --git a/experiments/standalone/Soil/water_conservation.jl b/experiments/standalone/Soil/water_conservation.jl index 595f2e77d2..c529e77674 100644 --- a/experiments/standalone/Soil/water_conservation.jl +++ b/experiments/standalone/Soil/water_conservation.jl @@ -11,6 +11,7 @@ using ClimaLand using ClimaLand.Soil using ClimaLand.Domains: Column import ClimaUtilities.OutputPathGenerator: generate_output_path +import ClimaLand.Simulations: LandSimulation, solve! rmse(v1, v2) = sqrt(mean((v1 .- v2) .^ 2)) @@ -81,43 +82,28 @@ for FT in (Float32, Float64) sources = sources, ) - exp_tendency! = make_exp_tendency(soil) - set_initial_cache! = make_set_initial_cache(soil) - imp_tendency! = make_imp_tendency(soil) - jacobian! = make_jacobian(soil) rmses = Array{FT}(undef, length(dts)) mass_errors = Array{FT}(undef, length(dts)) + function set_ic!(Y, p, t0, model) + @. Y.soil.ϑ_l = FT(0.24) + end for i in eachindex(dts) dt = dts[i] - - Y, p, coords = initialize(soil) - @. Y.soil.ϑ_l = FT(0.24) - set_initial_cache!(p, Y, FT(0.0)) + simulation = LandSimulation( + t_start, + t_end, + dt, + soil; + set_ic! = set_ic!, + timestepper = ode_algo, + solver_kwargs = (; saveat = collect(t_start:dt:t_end)), + ) + p = simulation._integrator.p p_init = deepcopy(p) mass_start = p_init.soil.total_water - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, - ) - - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t_start, t_end), - p, - ) + sol = solve!(simulation) - sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - saveat = collect(t_start:dt:t_end), - ) # Check that simulation still has correct float type @assert eltype(sol.u[end].soil) == FT @@ -215,44 +201,27 @@ for FT in (Float32, Float64) sources = sources, ) - exp_tendency! = make_exp_tendency(soil_dirichlet) - set_initial_cache! = make_set_initial_cache(soil_dirichlet) - imp_tendency! = make_imp_tendency(soil_dirichlet) - jacobian! = make_jacobian(soil_dirichlet) - update_aux! = make_update_aux(soil_dirichlet) - rmses_dirichlet = Array{FT}(undef, length(dts)) mass_errors_dirichlet = Array{FT}(undef, length(dts)) for i in eachindex(dts) dt = dts[i] - Y, p, coords = initialize(soil_dirichlet) - @. Y.soil.ϑ_l = FT(0.24) - set_initial_cache!(p, Y, FT(0.0)) + function set_ic!(Y, p, t0, model) + @. Y.soil.ϑ_l = FT(0.24) + end + simulation = LandSimulation( + t_start, + t_end, + dt, + soil_dirichlet; + set_ic! = set_ic!, + timestepper = ode_algo, + solver_kwargs = (; saveat = collect(t_start:dt:t_end)), + ) + p = simulation._integrator.p p_init = deepcopy(p) mass_start = p_init.soil.total_water - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, - ) - - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t_start, t_end), - p, - ) - sol = SciMLBase.solve( - prob, - ode_algo; - dt = dt, - adaptive = false, - saveat = Array(t_start:dt:t_end), - ) + sol = solve!(simulation) # Check that simulation still has correct float type @assert eltype(sol.u[end].soil) == FT diff --git a/experiments/standalone/Soil/water_energy_conservation.jl b/experiments/standalone/Soil/water_energy_conservation.jl index 8940e6f187..91976c5341 100644 --- a/experiments/standalone/Soil/water_energy_conservation.jl +++ b/experiments/standalone/Soil/water_energy_conservation.jl @@ -11,6 +11,7 @@ import ClimaParams as CP using ClimaLand using ClimaLand.Soil using ClimaLand.Domains: Column +import ClimaLand.Simulations: LandSimulation, solve! import ClimaUtilities.OutputPathGenerator: generate_output_path rmse(v1, v2) = sqrt(mean((v1 .- v2) .^ 2)) @@ -61,11 +62,6 @@ soil = Soil.EnergyHydrology{FT}(; sources = sources, ); -exp_tendency! = make_exp_tendency(soil); -imp_tendency! = make_imp_tendency(soil); -jacobian! = ClimaLand.make_jacobian(soil); - - function init_soil!(Y, z, Trange, params) ν = params.ν θ_r = params.θ_r @@ -101,7 +97,12 @@ function init_soil!(Y, z, Trange, params) params.earth_param_set, ) end -set_initial_cache! = make_set_initial_cache(soil); +function make_set_ic(z, Trange) + let z = z, Trange = Trange + (Y, p, t0, model) -> init_soil!(Y, z, Trange, model.parameters) + end +end + stepper = CTS.ARS111() err = (FT == Float64) ? 1e-8 : 1e-4 @@ -140,26 +141,17 @@ savedir = generate_output_path( for experiment in [no_phase_change, phase_change] (; t0, tf, dt_ref, Trange, dts, exp_name) = experiment - Y, p, coords = initialize(soil) - init_soil!(Y, coords.subsurface.z, Trange, soil.parameters) - set_initial_cache!(p, Y, t0) - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, + set_ic! = make_set_ic(soil.domain.fields.z, Trange) + ref_simulation = LandSimulation( + t0, + tf, + dt_ref, + soil; + set_ic! = set_ic!, + timestepper = ode_algo, + solver_kwargs = (; saveat = [t0, tf]), ) - - 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_sol = SciMLBase.solve(prob, ode_algo; dt = dt_ref, saveat = [t0, tf]) - + ref_sol = solve!(ref_simulation) # Now try several dts rmses_water = Array{FT}(undef, length(dts)) rmses_energy = Array{FT}(undef, length(dts)) @@ -168,29 +160,23 @@ for experiment in [no_phase_change, phase_change] for i in eachindex(dts) dt = dts[i] - Y, p, coords = initialize(soil) - init_soil!(Y, coords.subsurface.z, Trange, soil.parameters) - set_initial_cache!(p, Y, t0) + simulation = LandSimulation( + t0, + tf, + dt, + soil; + set_ic! = set_ic!, + timestepper = ode_algo, + solver_kwargs = (; saveat = [t0, tf]), + user_callbacks = (), + ) + p = simulation._integrator.p p_init = deepcopy(p) energy_start = Array(parent(p_init.soil.total_energy))[1] mass_start = Array(parent(p_init.soil.total_water))[1] - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, - ) - prob = SciMLBase.ODEProblem( - CTS.ClimaODEFunction( - T_exp! = exp_tendency!, - T_imp! = SciMLBase.ODEFunction(imp_tendency!; jac_kwargs...), - dss! = ClimaLand.dss!, - ), - Y, - (t0, tf), - p, - ) - sol = SciMLBase.solve(prob, ode_algo; dt = dt, saveat = [t0, tf]) + sol = solve!(simulation) # Calculate water mass balance over entire simulation mass_end = Array(parent(p.soil.total_water))[1] t_sim = sol.t[end] - sol.t[1] diff --git a/experiments/standalone/Vegetation/no_vegetation.jl b/experiments/standalone/Vegetation/no_vegetation.jl index 716cfcad13..a8ef1a8e40 100644 --- a/experiments/standalone/Vegetation/no_vegetation.jl +++ b/experiments/standalone/Vegetation/no_vegetation.jl @@ -17,6 +17,7 @@ using ClimaLand.Canopy.PlantHydraulics import ClimaLand import ClimaLand.Parameters as LP import ClimaUtilities.OutputPathGenerator: generate_output_path +import ClimaLand.Simulations: LandSimulation, solve! using DelimitedFiles import ClimaLand.FluxnetSimulations as FluxnetSimulations @@ -30,6 +31,9 @@ land_domain = Point(; z_sfc = FT(0.0), longlat = (long, lat)) atmos_h = FT(32) site_ID = "US-MOz" start_date = DateTime(2010) + Hour(time_offset) +N_days = 10 +end_date = start_date + Day(N_days) +dt = 225.0; # Prescribed forcing from Fluxnet data (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet( @@ -50,65 +54,49 @@ ground = PrescribedGroundConditions{FT}(; forcing = (; atmos, radiation, ground); function fakeLAIfunction(t) - if t < 30 * 24 * 3600 + if FT(t) < 30 * 24 * 3600 0.0 - elseif t < (364 - 30) * 24 * 3600.0 - max(2.0 * sin(2 * π / (730 * 24 * 3600) * (t - 30 * 24 * 3600)), 0) + elseif FT(t) < (364 - 30) * 24 * 3600.0 + max(2.0 * sin(2 * π / (730 * 24 * 3600) * (FT(t) - 30 * 24 * 3600)), 0) else 0.0 end end -LAI = TimeVaryingInput(fakeLAIfunction) -# Set up plant hydraulics with no vegetation +LAI = TimeVaryingInput(fakeLAIfunction) SAI = RAI = FT(0) hydraulics = PlantHydraulicsModel{FT}(land_domain, LAI; SAI, RAI); - -# Construct a CanopyModel with default parameters canopy = ClimaLand.Canopy.CanopyModel{FT}( land_domain, forcing, LAI, earth_param_set; hydraulics, -); +) + -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.FieldMatrixWithSolver(Y), Wfact = jacobian!); -ψ_leaf_0 = FT(-2e5 / 9800) -(; retention_model, ν, S_s) = canopy.hydraulics.parameters; -S_l_ini = inverse_water_retention_curve(retention_model, ψ_leaf_0, ν, S_s) -Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini) -t0 = 0.0 -N_days = 10 -tf = t0 + 3600 * 24 * N_days -dt = 225.0; -evaluate!(Y.canopy.energy.T, atmos.T, t0) -set_initial_cache! = make_set_initial_cache(canopy) -set_initial_cache!(p, Y, t0); +ψ_leaf_0 = FT(-2e5 / 9800) +(; retention_model, ν, S_s) = canopy.hydraulics.parameters; +S_l_ini = inverse_water_retention_curve(retention_model, ψ_leaf_0, ν, S_s) +function set_ic!(Y, p, t0, model) + Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini) + evaluate!(Y.canopy.energy.T, atmos.T, t0) +end n = 16 -saveat = Array(t0:(n * dt):tf) +saveat = Array(start_date:Second(n * dt):end_date) sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); -updateat = Array(t0:1800:tf) -drivers = ClimaLand.get_drivers(canopy) -updatefunc = ClimaLand.make_update_drivers(drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); +updateat = Array(start_date:Second(1800):end_date) # Set up timestepper timestepper = CTS.ARS111(); @@ -120,18 +108,18 @@ 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, -); - -sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); +simulation = LandSimulation( + start_date, + end_date, + dt, + canopy; + set_ic! = set_ic!, + user_callbacks = (saving_cb,), + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, +) +sol = solve!(simulation) savedir = generate_output_path("experiments/standalone/Vegetation/no_vegetation"); @@ -173,33 +161,34 @@ Tr = [ parent(sv.saveval[k].canopy.turbulent_fluxes.transpiration)[1] for k in 1:length(sol.t) ] +times = FT.(sol.t) ./ 24 ./ 3600 fig = Figure() ax = Axis(fig[1, 1], xlabel = "Time (days)", ylabel = "Temperature (K)") -lines!(ax, sol.t ./ 24 ./ 3600, T, label = "Canopy") -lines!(ax, sol.t ./ 24 ./ 3600, T_atmos, label = "Atmos") +lines!(ax, times, T, label = "Canopy") +lines!(ax, times, T_atmos, label = "Atmos") axislegend(ax) ax = Axis(fig[2, 1], xlabel = "Time (days)", ylabel = "Volumetric Water") -lines!(ax, sol.t ./ 24 ./ 3600, ϑ, label = "Canopy") +lines!(ax, times, ϑ, label = "Canopy") axislegend(ax) ax = Axis(fig[3, 1], xlabel = "Time (days)", ylabel = "LAI") -lines!(ax, sol.t ./ 24 ./ 3600, fakeLAIfunction.(sol.t), label = "Canopy") +lines!(ax, times, fakeLAIfunction.(sol.t), label = "Canopy") axislegend(ax) save(joinpath(savedir, "no_veg_state.png"), fig) fig2 = Figure() ax = Axis(fig2[1, 1], xlabel = "Time (days)", ylabel = "Energy Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, SW_n, label = "SW") -lines!(ax, sol.t ./ 24 ./ 3600, LW_n, label = "LW") -lines!(ax, sol.t ./ 24 ./ 3600, SHF, label = "SHF") -lines!(ax, sol.t ./ 24 ./ 3600, LHF, label = "LHF") -lines!(ax, sol.t ./ 24 ./ 3600, RE, label = "RE") +lines!(ax, times, SW_n, label = "SW") +lines!(ax, times, LW_n, label = "LW") +lines!(ax, times, SHF, label = "SHF") +lines!(ax, times, LHF, label = "LHF") +lines!(ax, times, RE, label = "RE") axislegend(ax) ax = Axis(fig2[2, 1], xlabel = "Time (days)", ylabel = "Water Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, Tr, label = "Transpiration") -lines!(ax, sol.t ./ 24 ./ 3600, R, label = "R") +lines!(ax, times, Tr, label = "Transpiration") +lines!(ax, times, R, label = "R") axislegend(ax) ax = Axis(fig2[3, 1], xlabel = "Time (days)", ylabel = "Carbon Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, GPP, label = "GPP") -lines!(ax, sol.t ./ 24 ./ 3600, resp, label = "Respiration") +lines!(ax, times, GPP, label = "GPP") +lines!(ax, times, resp, label = "Respiration") axislegend(ax) save(joinpath(savedir, "no_veg_fluxes.png"), fig2) diff --git a/experiments/standalone/Vegetation/timestep_test.jl b/experiments/standalone/Vegetation/timestep_test.jl index 40794cd9a0..ccd9976cb2 100644 --- a/experiments/standalone/Vegetation/timestep_test.jl +++ b/experiments/standalone/Vegetation/timestep_test.jl @@ -55,6 +55,7 @@ using ClimaLand.Domains: Point using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics import ClimaLand +import ClimaLand.Simulations: LandSimulation, solve! import ClimaLand.Parameters as LP using DelimitedFiles import ClimaLand.FluxnetSimulations as FluxnetSimulations @@ -69,11 +70,9 @@ long = FT(-92.2000) # degree land_domain = Point(; z_sfc = FT(0.0), longlat = (long, lat)) atmos_h = FT(32) site_ID = "US-MOz" -start_date = DateTime(2010) + Hour(time_offset) # UTC -seconds_per_day = 3600 * 24.0 -t0 = 150seconds_per_day +start_date = DateTime(2010) + Hour(time_offset) + Day(150) N_days = 20.0 -tf = t0 + N_days * seconds_per_day + 80 +end_date = start_date + Day(N_days) + Second(80) # Get prescribed atmospheric and radiation forcing (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet( @@ -97,8 +96,8 @@ forcing = (; atmos, radiation, ground); surface_space = land_domain.space.surface modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(; context = ClimaComms.context(surface_space), - start_date = start_date + Second(t0), - end_date = start_date + Second(t0) + Second(tf), + start_date = start_date, + end_date = end_date, ) LAI = ClimaLand.prescribed_lai_modis( modis_lai_ncdata_path, @@ -118,12 +117,6 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}( energy, ); -exp_tendency! = make_exp_tendency(canopy) -imp_tendency! = make_imp_tendency(canopy) -jacobian! = make_jacobian(canopy) -drivers = ClimaLand.get_drivers(canopy) -updatefunc = ClimaLand.make_update_drivers(drivers) - (; retention_model, ν, S_s) = canopy.hydraulics.parameters; ψ_leaf_0 = FT(-2e5 / 9800) ψ_stem_0 = FT(-1e5 / 9800) @@ -135,8 +128,6 @@ S_l_ini = S_s, ) -set_initial_cache! = make_set_initial_cache(canopy) - timestepper = CTS.ARS111(); ode_algo = CTS.IMEXAlgorithm( timestepper, @@ -145,7 +136,10 @@ ode_algo = CTS.IMEXAlgorithm( update_j = CTS.UpdateEvery(CTS.NewNewtonIteration), ), ); - +function set_ic!(Y, p, t0, model) + Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1]) + evaluate!(Y.canopy.energy.T, atmos.T, t0) +end ref_dt = 6.0 dts = [ref_dt, 12.0, 48.0, 225.0, 450.0, 900.0, 1800.0, 3600.0] @@ -160,37 +154,21 @@ times = [] for dt in dts @info dt - # Initialize model before each simulation - Y, p, coords = ClimaLand.initialize(canopy) - jac_kwargs = (; - jac_prototype = ClimaLand.FieldMatrixWithSolver(Y), - Wfact = jacobian!, + # Create update times for driver and saving callback + saveat = vcat(Array(start_date:Second(3 * 3600):end_date), end_date) + updateat = vcat(Array(start_date:Second(3 * 3600):end_date), end_date) + simulation = LandSimulation( + start_date, + end_date, + dt, + canopy; + set_ic! = set_ic!, + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, + user_callbacks = (), ) - - # Set initial conditions - Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini[1]) - evaluate!(Y.canopy.energy.T, atmos.T, t0) - set_initial_cache!(p, Y, t0) - - # Create callback for driver updates - saveat = vcat(Array(t0:(3 * 3600):tf), tf) - updateat = vcat(Array(t0:(3 * 3600):tf), tf) - cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) - - # Solve simulation - 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) - + @time sol = solve!(simulation) # Save results for comparison if dt == ref_dt global ref_T = @@ -204,7 +182,7 @@ for dt in dts push!(p99_err, FT(percentile(ΔT, 99))) push!(sol_err, ΔT[end]) push!(T_states, T) - push!(times, sol.t) + push!(times, float.(sol.t)) end end savedir = generate_output_path( @@ -212,7 +190,7 @@ savedir = generate_output_path( ); # Create plot with statistics -sim_time = round((tf - t0) / 3600, digits = 2) # simulation length in hours +sim_time = round(Dates.value(Second(end_date - start_date)) / 3600, digits = 2) # simulation length in hours # Compare T state computed with small vs large dt fig = Figure() diff --git a/experiments/standalone/Vegetation/varying_lai.jl b/experiments/standalone/Vegetation/varying_lai.jl index bba11617ee..6e510cac1b 100644 --- a/experiments/standalone/Vegetation/varying_lai.jl +++ b/experiments/standalone/Vegetation/varying_lai.jl @@ -16,6 +16,7 @@ using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics import ClimaLand import ClimaLand.Parameters as LP +import ClimaLand.Simulations: LandSimulation, solve! using DelimitedFiles import ClimaLand.FluxnetSimulations as FluxnetSimulations import ClimaUtilities.OutputPathGenerator: generate_output_path @@ -29,6 +30,9 @@ land_domain = Point(; z_sfc = FT(0.0), longlat = (long, lat)) atmos_h = FT(32) site_ID = "US-MOz" start_date = DateTime(2010) + Hour(time_offset) +N_days = 364 +end_date = start_date + Day(N_days) +dt = 225.0; (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet( site_ID, lat, @@ -47,10 +51,10 @@ ground = PrescribedGroundConditions{FT}(; forcing = (; atmos, radiation, ground); function fakeLAIfunction(t) - if t < 30 * 24 * 3600 + if FT(t) < 30 * 24 * 3600 0.0 - elseif t < (364 - 30) * 24 * 3600.0 - max(2.0 * sin(2 * π / (730 * 24 * 3600) * (t - 30 * 24 * 3600)), 0) + elseif FT(t) < (364 - 30) * 24 * 3600.0 + max(2.0 * sin(2 * π / (730 * 24 * 3600) * (FT(t) - 30 * 24 * 3600)), 0) else 0.0 end @@ -73,41 +77,26 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}( hydraulics, ); -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.FieldMatrixWithSolver(Y), Wfact = jacobian!); - (; retention_model, ν, S_s) = canopy.hydraulics.parameters; ψ_leaf_0 = FT(-2e5 / 9800) S_l_ini = inverse_water_retention_curve(retention_model, ψ_leaf_0, ν, S_s) -Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini) -t0 = 0.0 -N_days = 364 -tf = t0 + 3600 * 24 * N_days -dt = 225.0; -evaluate!(Y.canopy.energy.T, atmos.T, t0) -set_initial_cache! = make_set_initial_cache(canopy) -set_initial_cache!(p, Y, t0); +function set_ic!(Y, p, t0, model) + Y.canopy.hydraulics.ϑ_l.:1 .= augmented_liquid_fraction.(ν, S_l_ini) + evaluate!(Y.canopy.energy.T, atmos.T, t0) +end n = 16 -saveat = Array(t0:(n * dt):tf) +saveat = Array(start_date:Second(n * dt):end_date) sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); - -updateat = Array(t0:1800:tf) +updateat = Array(start_date:Second(1800):end_date) drivers = ClimaLand.get_drivers(canopy) -updatefunc = ClimaLand.make_update_drivers(drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); # Set up timestepper timestepper = CTS.ARS111(); @@ -119,18 +108,18 @@ 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, -); - -sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); +simulation = LandSimulation( + start_date, + end_date, + dt, + canopy; + set_ic! = set_ic!, + user_callbacks = (saving_cb,), + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, +) +sol = solve!(simulation) savedir = generate_output_path("experiments/standalone/Vegetation/varying_lai"); T = [parent(sol.u[k].canopy.energy.T)[1] for k in 1:length(sol.t)] @@ -171,33 +160,34 @@ Tr = [ parent(sv.saveval[k].canopy.turbulent_fluxes.transpiration)[1] for k in 1:length(sol.t) ] +times = float.(sol.t) ./ 24 ./ 3600 fig = Figure() ax = Axis(fig[1, 1], xlabel = "Time (days)", ylabel = "Temperature (K)") -lines!(ax, sol.t ./ 24 ./ 3600, T, label = "Canopy") -lines!(ax, sol.t ./ 24 ./ 3600, T_atmos, label = "Atmos") +lines!(ax, times, T, label = "Canopy") +lines!(ax, times, T_atmos, label = "Atmos") axislegend(ax) ax = Axis(fig[2, 1], xlabel = "Time (days)", ylabel = "Volumetric Water") -lines!(ax, sol.t ./ 24 ./ 3600, ϑ, label = "Canopy") +lines!(ax, times, ϑ, label = "Canopy") axislegend(ax) ax = Axis(fig[3, 1], xlabel = "Time (days)", ylabel = "LAI") -lines!(ax, sol.t ./ 24 ./ 3600, fakeLAIfunction.(sol.t), label = "Canopy") +lines!(ax, times, fakeLAIfunction.(sol.t), label = "Canopy") axislegend(ax) save(joinpath(savedir, "varying_lai_state.png"), fig) fig2 = Figure() ax = Axis(fig2[1, 1], xlabel = "Time (days)", ylabel = "Energy Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, SW_n, label = "SW_n") -lines!(ax, sol.t ./ 24 ./ 3600, LW_n, label = "LW_n") -lines!(ax, sol.t ./ 24 ./ 3600, SHF, label = "SHF") -lines!(ax, sol.t ./ 24 ./ 3600, LHF, label = "LHF") -lines!(ax, sol.t ./ 24 ./ 3600, RE, label = "RE") +lines!(ax, times, SW_n, label = "SW_n") +lines!(ax, times, LW_n, label = "LW_n") +lines!(ax, times, SHF, label = "SHF") +lines!(ax, times, LHF, label = "LHF") +lines!(ax, times, RE, label = "RE") axislegend(ax) ax = Axis(fig2[2, 1], xlabel = "Time (days)", ylabel = "Water Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, Tr, label = "Transpiration") -lines!(ax, sol.t ./ 24 ./ 3600, R, label = "R") +lines!(ax, times, Tr, label = "Transpiration") +lines!(ax, times, R, label = "R") axislegend(ax) ax = Axis(fig2[3, 1], xlabel = "Time (days)", ylabel = "Carbon Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, GPP, label = "GPP") -lines!(ax, sol.t ./ 24 ./ 3600, resp, label = "Respiration") +lines!(ax, times, GPP, label = "GPP") +lines!(ax, times, resp, label = "Respiration") axislegend(ax) save(joinpath(savedir, "varying_lai_fluxes.png"), fig2) diff --git a/experiments/standalone/Vegetation/varying_lai_with_stem.jl b/experiments/standalone/Vegetation/varying_lai_with_stem.jl index 3f859b7aa8..76f49d23f3 100644 --- a/experiments/standalone/Vegetation/varying_lai_with_stem.jl +++ b/experiments/standalone/Vegetation/varying_lai_with_stem.jl @@ -16,6 +16,7 @@ using ClimaLand.Canopy using ClimaLand.Canopy.PlantHydraulics import ClimaLand import ClimaLand.Parameters as LP +import ClimaLand.Simulations: LandSimulation, solve! import ClimaUtilities.OutputPathGenerator: generate_output_path using DelimitedFiles import ClimaLand.FluxnetSimulations as FluxnetSimulations @@ -29,6 +30,9 @@ land_domain = Point(; z_sfc = FT(0.0), longlat = (long, lat)) atmos_h = FT(32) site_ID = "US-MOz" start_date = DateTime(2010) + Hour(time_offset) +N_days = 60 +end_date = start_date + Day(N_days) +dt = FT(225); (; atmos, radiation) = FluxnetSimulations.prescribed_forcing_fluxnet( site_ID, lat, @@ -39,6 +43,7 @@ start_date = DateTime(2010) + Hour(time_offset) earth_param_set, FT, ) + ground = PrescribedGroundConditions{FT}(; α_PAR = FT(0.2), α_NIR = FT(0.4), @@ -47,10 +52,10 @@ ground = PrescribedGroundConditions{FT}(; forcing = (; atmos, radiation, ground); function fakeLAIfunction2(t) - if t < 10 * 24 * 3600 + if FT(t) < 10 * 24 * 3600 0.0 - elseif t < (60 - 10) * 24 * 3600.0 - max(2.0 * sin(2 * π / (730 * 24 * 3600) * (t - 10 * 24 * 3600)), 0) + elseif FT(t) < (60 - 10) * 24 * 3600.0 + max(2.0 * sin(2 * π / (730 * 24 * 3600) * (FT(t) - 10 * 24 * 3600)), 0) else 0.0 end @@ -105,13 +110,6 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}( ); -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.FieldMatrixWithSolver(Y), Wfact = jacobian!); - ψ_leaf_0 = FT(-2e5 / 9800) ψ_stem_0 = FT(-1e5 / 9800) @@ -123,31 +121,22 @@ 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]) - -t0 = 0.0 -N_days = 60 -tf = t0 + 3600 * 24 * N_days -dt = 225.0; -evaluate!(Y.canopy.energy.T, atmos.T, t0) -set_initial_cache! = make_set_initial_cache(canopy) -set_initial_cache!(p, Y, t0); +function set_ic!(Y, p, t0, model) + 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) +end n = 16 -saveat = Array(t0:(n * dt):tf) +saveat = Array(start_date:Second(n * dt):end_date) sv = (; - t = Array{Float64}(undef, length(saveat)), + t = Array{DateTime}(undef, length(saveat)), saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat); -updateat = Array(t0:1800:tf) -drivers = ClimaLand.get_drivers(canopy) -updatefunc = ClimaLand.make_update_drivers(drivers) -driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) -cb = SciMLBase.CallbackSet(driver_cb, saving_cb); +updateat = Array(start_date:Second(1800):end_date) # Set up timestepper timestepper = CTS.ARS111(); @@ -159,19 +148,18 @@ 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, -); - -sol = SciMLBase.solve(prob, ode_algo; dt = dt, callback = cb, saveat = saveat); - +simulation = LandSimulation( + start_date, + end_date, + dt, + canopy; + set_ic! = set_ic!, + user_callbacks = (saving_cb,), + updateat = updateat, + solver_kwargs = (; saveat = deepcopy(saveat)), + timestepper = ode_algo, +) +sol = solve!(simulation) savedir = generate_output_path( "experiments/standalone/Vegetation/varying_lai_with_stem", ); @@ -216,36 +204,37 @@ Tr = [ parent(sv.saveval[k].canopy.turbulent_fluxes.transpiration)[1] for k in 1:length(sol.t) ] +times = float.(sol.t) ./ 24 ./ 3600 fig = Figure() ax = Axis(fig[1, 1], xlabel = "Time (days)", ylabel = "Temperature (K)") -lines!(ax, sol.t ./ 24 ./ 3600, T, label = "Canopy") -lines!(ax, sol.t ./ 24 ./ 3600, T_atmos, label = "Atmos") +lines!(ax, times, T, label = "Canopy") +lines!(ax, times, T_atmos, label = "Atmos") axislegend(ax) ax = Axis(fig[2, 1], xlabel = "Time (days)", ylabel = "Volumetric Water") -lines!(ax, sol.t ./ 24 ./ 3600, ϑ_l, label = "Leaf") -lines!(ax, sol.t ./ 24 ./ 3600, ϑ_s, label = "Stem") +lines!(ax, times, ϑ_l, label = "Leaf") +lines!(ax, times, ϑ_s, label = "Stem") axislegend(ax) ax = Axis(fig[3, 1], xlabel = "Time (days)", ylabel = "LAI") -lines!(ax, sol.t ./ 24 ./ 3600, fakeLAIfunction2.(sol.t), label = "Canopy") +lines!(ax, times, fakeLAIfunction2.(sol.t), label = "Canopy") axislegend(ax) save(joinpath(savedir, "varying_lai_with_stem_state.png"), fig) fig2 = Figure() ax = Axis(fig2[1, 1], xlabel = "Time (days)", ylabel = "Energy Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, SW_n, label = "SW_n") -lines!(ax, sol.t ./ 24 ./ 3600, LW_n, label = "LW_n") -lines!(ax, sol.t ./ 24 ./ 3600, SHF, label = "SHF") -lines!(ax, sol.t ./ 24 ./ 3600, LHF, label = "LHF") -lines!(ax, sol.t ./ 24 ./ 3600, RE, label = "RE") +lines!(ax, times, SW_n, label = "SW_n") +lines!(ax, times, LW_n, label = "LW_n") +lines!(ax, times, SHF, label = "SHF") +lines!(ax, times, LHF, label = "LHF") +lines!(ax, times, RE, label = "RE") axislegend(ax) ax = Axis(fig2[2, 1], xlabel = "Time (days)", ylabel = "Water Fluxes") ylims!(ax, (0, 1e-7)) -lines!(ax, sol.t ./ 24 ./ 3600, Tr, label = "Transpiration") -lines!(ax, sol.t ./ 24 ./ 3600, R, label = "R") -lines!(ax, sol.t ./ 24 ./ 3600, R_stem_leaf, label = "R_stem_leaf") +lines!(ax, times, Tr, label = "Transpiration") +lines!(ax, times, R, label = "R") +lines!(ax, times, R_stem_leaf, label = "R_stem_leaf") axislegend(ax) ax = Axis(fig2[3, 1], xlabel = "Time (days)", ylabel = "Carbon Fluxes") -lines!(ax, sol.t ./ 24 ./ 3600, GPP, label = "GPP") -lines!(ax, sol.t ./ 24 ./ 3600, resp, label = "Respiration") +lines!(ax, times, GPP, label = "GPP") +lines!(ax, times, resp, label = "Respiration") axislegend(ax) save(joinpath(savedir, "varying_lai_with_stem_fluxes.png"), fig2) diff --git a/ext/LandSimulationVisualizationExt.jl b/ext/LandSimulationVisualizationExt.jl index 970da0be80..352bb04c7d 100644 --- a/ext/LandSimulationVisualizationExt.jl +++ b/ext/LandSimulationVisualizationExt.jl @@ -6,6 +6,7 @@ import ClimaAnalysis import ClimaAnalysis.Visualize as viz using CairoMakie import GeoMakie +import ClimaUtilities.TimeManager: ITime using Dates import NCDatasets using ClimaLand diff --git a/ext/fluxnet_simulations/forcing.jl b/ext/fluxnet_simulations/forcing.jl index 50959cb942..8789068f2e 100644 --- a/ext/fluxnet_simulations/forcing.jl +++ b/ext/fluxnet_simulations/forcing.jl @@ -247,21 +247,46 @@ function FluxnetSimulations.get_data_dt(site_ID) end """ - get_data_dates(site_ID, hour_offset_from_UTC) + get_data_dates( + site_ID, + hour_offset_from_UTC; + duration::Union{Nothing, Period} = nothing, + start_offset::Period = Second(0), + ) A helper function to get the first and last dates, in UTC, for which we have Fluxnet data at `site_ID`, given the offset in hours of local time -from UTC. +from UTC. If `duration` is provided, it is used to determine the end date, +otherwise the end date is the last date in the data. The `start_offset` is +added to the start date, and must be non-negative. """ -function FluxnetSimulations.get_data_dates(site_ID, hour_offset_from_UTC) +function FluxnetSimulations.get_data_dates( + site_ID, + hour_offset_from_UTC; + duration::Union{Nothing, Period} = nothing, + start_offset::Period = Second(0), +) fluxnet_csv_path = ClimaLand.Artifacts.experiment_fluxnet_data_path(site_ID) (data, columns) = readdlm(fluxnet_csv_path, ','; header = true) local_datetime = DateTime.(string.(Int.(data[:, 1])), "yyyymmddHHMM") UTC_datetime = local_datetime .+ Dates.Hour(hour_offset_from_UTC) - return extrema(UTC_datetime) + earliest_date, latest_date = extrema(UTC_datetime) + Dates.value(start_offset) < 0 && error("start_offset must be non-negative") + if !isnothing(duration) && Dates.value(duration) < 0 + error("If duration is not provided, it must be non-negative.") + end + duration_available_ms = latest_date - earliest_date + start_date = earliest_date + start_offset + end_date = isnothing(duration) ? latest_date : start_date + duration + if !isnothing(duration) && (end_date > latest_date) + error( + "The sum of the requested duration of $duration and start_offset of $start_offset \ + is greater than the available $duration_available_ms of data.", + ) + end + return (start_date, end_date) end - """ snow_precip_fraction(T_in_C, VPD_in_hPa; thermo_params) diff --git a/ext/land_sim_vis/plotting_utils.jl b/ext/land_sim_vis/plotting_utils.jl index e67726494a..befa230b56 100644 --- a/ext/land_sim_vis/plotting_utils.jl +++ b/ext/land_sim_vis/plotting_utils.jl @@ -369,7 +369,7 @@ function LandSimVis.make_diurnal_timeseries( dn, ) save_Δt = model_time[2] - model_time[1] # in seconds - model_dates = Second.(model_time) .+ start_date + model_dates = Second.(float.(model_time)) .+ start_date spinup_idx = findfirst(spinup_date .<= model_dates) hour_of_day, model_diurnal_cycle = compute_diurnal_cycle( model_dates[spinup_idx:end], @@ -486,7 +486,7 @@ function LandSimVis.make_timeseries( dn, ) save_Δt = model_time[2] - model_time[1] # in seconds - model_dates = Second.(model_time) .+ start_date + model_dates = Second.(float.(model_time)) .+ start_date spinup_idx = findfirst(spinup_date .<= model_dates) hour_of_day, model_diurnal_cycle = compute_diurnal_cycle( model_dates[spinup_idx:end], diff --git a/src/Artifacts.jl b/src/Artifacts.jl index 60041a5e5c..6a9053ce30 100644 --- a/src/Artifacts.jl +++ b/src/Artifacts.jl @@ -83,7 +83,7 @@ function era5_land_forcing_data2008_path(; context = nothing, lowres = false) return hires_path catch @warn( - "High resolution ERA5 forcing not available locally; downloading and using low resolution data instead." + "High resolution ERA5 forcing not available locally; using low resolution data instead." ) return lowres_path end diff --git a/src/diagnostics/Diagnostics.jl b/src/diagnostics/Diagnostics.jl index c62c8f807f..d1ce2e93d1 100644 --- a/src/diagnostics/Diagnostics.jl +++ b/src/diagnostics/Diagnostics.jl @@ -18,7 +18,7 @@ import ..Soil: EnergyHydrology import ..Canopy: medlyn_conductance, medlyn_term, moisture_stress import ..Domains: top_center_to_surface, AbstractDomain, SphericalShell, HybridBox - +import ClimaLand import ..heaviside import ClimaDiagnostics: @@ -32,6 +32,7 @@ import ClimaDiagnostics.Writers: HDF5Writer, NetCDFWriter, DictWriter import ClimaCore.Fields: zeros, field_values import ClimaCore.Operators: column_integral_definite! +import ClimaUtilities.TimeManager: ITime, date export close_output_writers diff --git a/src/diagnostics/default_diagnostics.jl b/src/diagnostics/default_diagnostics.jl index 096b7e603e..1bfe032734 100644 --- a/src/diagnostics/default_diagnostics.jl +++ b/src/diagnostics/default_diagnostics.jl @@ -49,6 +49,23 @@ end include("standard_diagnostic_frequencies.jl") +# The default diagnostics currently require a start date because they use Dates.Period. +default_diagnostics( + model::ClimaLand.AbstractModel, + start_date::Union{Nothing, ITime{<:Any, <:Any, Nothing}}, + outdir, +) = [] + +function default_diagnostics(model::ClimaLand.AbstractModel, start_date, outdir) + # a starting date is required for default diagnostics + domain = ClimaLand.get_domain(model) + default_diagnostic_domain = + haskey(domain.space, :subsurface) ? domain.space.subsurface : + domain.space.surface + output_writer = NetCDFWriter(default_diagnostic_domain, outdir; start_date) + return default_diagnostics(model, start_date; output_writer) +end + # Bucket function default_diagnostics( land_model::BucketModel{FT}, @@ -401,3 +418,11 @@ function default_diagnostics( ) end end + +# fallback method if no default diagnostics are defined +default_diagnostics( + land_model::ClimaLand.AbstractModel, + start_date; + output_writer, + average_period = nothing, +) = [] diff --git a/src/shared_utilities/Domains.jl b/src/shared_utilities/Domains.jl index 2e87422406..94555443af 100644 --- a/src/shared_utilities/Domains.jl +++ b/src/shared_utilities/Domains.jl @@ -1213,7 +1213,16 @@ function landsea_mask( return binary_mask end -function landsea_mask(domain::Domains.AbstractDomain; kwargs...) +# Points and Columns do not have a horizontal dim, so a horizontal mask cannot be applied +landsea_mask(domain::Union{Point, Column}; kwargs...) = nothing + +function landsea_mask( + domain::Union{SphericalShell, SphericalSurface, HybridBox, Plane}; + kwargs..., +) + # HybridBox and Plane domains might not have longlat, which is needed for the mask + (hasproperty(domain, :longlat) && isnothing(domain.longlat)) && + return nothing # average_horizontal_resolution_degrees returns a tuple with the resolution # along the two directions, so we take the minimum resolution_degrees = minimum(average_horizontal_resolution_degrees(domain)) diff --git a/src/shared_utilities/drivers.jl b/src/shared_utilities/drivers.jl index a515bb77cf..89a762f65d 100644 --- a/src/shared_utilities/drivers.jl +++ b/src/shared_utilities/drivers.jl @@ -285,8 +285,12 @@ function default_zenith_angle( insol_params, ) where {T, LT} FT = eltype(latitude) - current_datetime = - T <: ITime ? date(t) : start_date + Dates.Second(round(t)) + current_datetime = if T <: ITime + # ITime may not have an epoch, so use start_date as fallback + isnothing(t.epoch) ? start_date + t.counter * t.period : date(t) + else + start_date + Dates.Second(round(t)) + end d, δ, η_UTC = FT.( Insolation.helper_instantaneous_zenith_angle( @@ -2060,7 +2064,12 @@ function empirical_diffuse_fraction( thermo_params, ) where {FT} if t isa ITime - DOY = Dates.dayofyear(date(t)) + # If t is an ITime, and it has an epoch, use that instead of start_date to compute + # the current simulation date + DOY = Dates.dayofyear( + isnothing(t.epoch) ? + start_date + Dates.Second(floor(Int64, float(t))) : date(t), + ) else DOY = Dates.dayofyear(start_date + Dates.Second(floor(Int64, t))) end diff --git a/src/shared_utilities/utils.jl b/src/shared_utilities/utils.jl index 3f89d6a1c7..2f9d9f261a 100644 --- a/src/shared_utilities/utils.jl +++ b/src/shared_utilities/utils.jl @@ -1,6 +1,6 @@ import ClimaCore import SciMLBase -import ClimaDiagnostics.Schedules: EveryCalendarDtSchedule +import ClimaDiagnostics.Schedules: EveryCalendarDtSchedule, EveryDtSchedule import Dates using ClimaComms @@ -9,6 +9,7 @@ import Interpolations import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput import ClimaUtilities.Regridders: InterpolationsRegridder import ClimaUtilities.OnlineLogging: WallTimeInfo, report_walltime +import ClimaUtilities.TimeManager: ITime, date export FTfromY, call_count_nans_state @@ -397,11 +398,12 @@ end This function is used by `NonInterpSavingCallback` to perform the saving. """ function (affect!::SavingAffect)(integrator) - while !isempty(affect!.saveat) && first(affect!.saveat) <= integrator.t + integrator_t = convert_t_saveat(affect!.saveat, integrator.t) + while !isempty(affect!.saveat) && first(affect!.saveat) <= integrator_t affect!.saveiter += 1 curr_t = popfirst!(affect!.saveat) # @assert curr_t == integrator.t - if curr_t == integrator.t + if curr_t == integrator_t affect!.saved_values.t[affect!.saveiter] = curr_t affect!.saved_values.saveval[affect!.saveiter] = deepcopy(integrator.p) @@ -409,6 +411,9 @@ function (affect!::SavingAffect)(integrator) end end +convert_t_saveat(saveat, t) = eltype(saveat)(t) +convert_t_saveat(saveat::Vector{T}, t::T) where {T} = t + """ saving_initialize(cb, u, t, integrator) @@ -418,7 +423,8 @@ initial values, don't pass the `initialize` argument to the `DiscreteCallback` constructor. """ function saving_initialize(cb, u, t, integrator) - (t in cb.affect!.saveat) && cb.affect!(integrator) + saveat = cb.affect!.saveat + (convert_t_saveat(saveat, t) in saveat) && cb.affect!(integrator) end """ @@ -457,7 +463,7 @@ This function returns a function with the type signature expected by called in the callback. This implementation simply checks if the current time is contained in the list of affect times used for the callback. """ -condition(saveat) = (_, t, _) -> t in saveat +condition(saveat) = (_, t, _) -> convert_t_saveat(saveat, t) in saveat function FTfromY(Y::ClimaCore.Fields.FieldVector) return eltype(Y) @@ -505,6 +511,10 @@ function isdivisible(dt_large::Dates.OtherPeriod, dt_small::Dates.OtherPeriod) return isinteger(dt_large / dt_small) end +function isdivisible(dt_large::ITime, dt_small::ITime) + return isinteger(dt_large / dt_small) +end + function isdivisible( dt_large::Union{Dates.Month, Dates.Year}, dt_small::Dates.FixedPeriod, @@ -559,7 +569,10 @@ end function call_count_nans_state( state::ClimaCore.Fields.Field, - space::ClimaCore.Spaces.AbstractSpectralElementSpace; + space::Union{ + ClimaCore.Spaces.AbstractSpectralElementSpace, + ClimaCore.Spaces.AbstractPointSpace, + }; mask = nothing, verbose = false, ) @@ -568,7 +581,10 @@ end function call_count_nans_state( state::ClimaCore.Fields.Field, - space::ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace; + space::Union{ + ClimaCore.Spaces.ExtrudedFiniteDifferenceSpace, + ClimaCore.Spaces.FiniteDifferenceSpace, + }; mask = nothing, verbose = false, ) @@ -628,7 +644,7 @@ initialized with the `start_date` to ensure that it is first called at the correct time. """ function NaNCheckCallback( - nancheck_frequency::Union{AbstractFloat, Dates.Period}, + nancheck_frequency::Union{AbstractFloat, Dates.Period, ITime}, start_date, dt; mask = nothing, @@ -652,7 +668,9 @@ function NaNCheckCallback( ) if !isnothing(dt) - dt_period = Dates.Millisecond(1000float(dt)) + dt_period = + typeof(dt) <: typeof(nancheck_frequency_period) ? dt : + Dates.Millisecond(1000float(dt)) if !isdivisible(nancheck_frequency_period, dt_period) @warn "Callback frequency ($(nancheck_frequency_period)) is not an integer multiple of dt $(dt_period)" end @@ -666,6 +684,27 @@ function NaNCheckCallback( SciMLBase.DiscreteCallback(cond, affect!) end +# no date version +function NaNCheckCallback( + nancheck_frequency::Union{AbstractFloat, ITime}, + start_date::Union{Nothing, ITime{<:Any, <:Any, Nothing}}, + dt; + mask = nothing, +) + schedule = EveryDtSchedule(nancheck_frequency) + if !isnothing(dt) + if !isdivisible(nancheck_frequency, dt) + @warn "Callback frequency ($(nancheck_frequency)) is not an integer multiple of dt $(dt)" + end + end + + cond = let schedule = schedule + (u, t, integrator) -> schedule(integrator) + end + affect! = (integrator) -> call_count_nans_state(integrator.u; mask) + SciMLBase.DiscreteCallback(cond, affect!) +end + """ ReportCallback(every_n_steps) diff --git a/src/simulations/Simulations.jl b/src/simulations/Simulations.jl index 53ea42a122..9552d38892 100644 --- a/src/simulations/Simulations.jl +++ b/src/simulations/Simulations.jl @@ -7,6 +7,7 @@ using Dates import ClimaUtilities.TimeManager: ITime, date import ClimaDiagnostics using ClimaLand +import ClimaLand.Domains: SphericalShell, HybridBox, SphericalSurface export step!, solve!, LandSimulation include("initial_conditions.jl") @@ -23,7 +24,7 @@ import ..Diagnostics: close_output_writers I <: SciMLBase.DEIntegrator, } -the ClimaLand LandSimulation struct, which specifies +the ClimaLand LandSimulation struct, which specifies - the discrete set of equations to solve (defined by the `model`); - the timestepping algorithm; - user callbacks (passed as a tuple) to be executed at specific times in the simulations; @@ -33,13 +34,13 @@ User callbacks are optional: examples currently include callbacks that estimate to solution and SYPD of the simulation as it runs, checkpoint the state, or check the solution for NaNs. Others can be added here. -Diagnostics are implemented as callbacks, and are also optional. -However, a default is provided. `diagnostics` is expected to be a +Diagnostics are implemented as callbacks, and are also optional. +However, a default is provided. `diagnostics` is expected to be a list of `ClimaDiagnostics.ScheduledDiagnostics`. Finally, the private field _required_callbacks consists of callbacks that are required for the simulation to run correctly. Currently, this includes the callbacks which update the atmospheric -forcing and update the LAI using prescribed data. +forcing and update the LAI using prescribed data. """ struct LandSimulation{ M <: ClimaLand.AbstractModel, @@ -59,11 +60,57 @@ struct LandSimulation{ _integrator::I end +""" + LandSimulation( + t0::ITime, + tf::ITime, + Δt::ITime, + model; + outdir = ".", + set_ic! = make_set_initial_state_from_file( + ClimaLand.Artifacts.soil_ic_2008_50m_path(; + context = ClimaComms.context(model), + ), + model, + ), + timestepper = ClimaTimeSteppers.IMEXAlgorithm( + ClimaTimeSteppers.ARS111(), + ClimaTimeSteppers.NewtonsMethod( + max_iters = 3, + update_j = ClimaTimeSteppers.UpdateEvery( + ClimaTimeSteppers.NewNewtonIteration, + ), + ), + ), + user_callbacks = ( + ClimaLand.NaNCheckCallback( + isnothing(t0.epoch) ? Δt * 10000 : Dates.Month(1), + t0, + Δt; + mask = ClimaLand.Domains.landsea_mask(ClimaLand.get_domain(model)), + ), + ClimaLand.ReportCallback(1000), + ), + diagnostics = ClimaLand.default_diagnostics(model, t0, outdir), + updateat = [promote(t0:(ITime(3600 * 3)):tf...)...], + solver_kwargs = (;), + ) + +Creates a `LandSimulation` object for `model` with the `_integrator` field initialized +with a problem from `t0` to `tf`, a time step of `Δt`, and with `timestepper` as the +time-stepping algorithm. During the construction process, the initial conditions +are set using `set_ic!`, the initial cache is set using `make_set_initial_cache(model)`, +and driver and diagnostic callbacks are created. If `updateat` is a vector, the drivers +will be updated at those times; if it is a single time, the drivers will be updated +at that interval. The order the callbacks are called while stepping the simulation is: +1. user_callbacks +2. driver update callback +3. diagnostics callback +""" function LandSimulation( - FT, - start_date::Dates.DateTime, - stop_date::Dates.DateTime, - Δt::AbstractFloat, + t0::ITime, + tf::ITime, + Δt::ITime, model; outdir = ".", set_ic! = make_set_initial_state_from_file( @@ -83,34 +130,27 @@ function LandSimulation( ), user_callbacks = ( ClimaLand.NaNCheckCallback( - Dates.Month(1), - start_date, - ITime(Δt, epoch = start_date), + isnothing(t0.epoch) ? Δt * 10000 : Dates.Month(1), + t0, + Δt; mask = ClimaLand.Domains.landsea_mask(ClimaLand.get_domain(model)), ), ClimaLand.ReportCallback(1000), ), - diagnostics = ClimaLand.default_diagnostics( - model, - start_date; - output_writer = ClimaDiagnostics.Writers.NetCDFWriter( - ClimaLand.get_domain(model).space.subsurface, - outdir; - start_date, - ), - ), + diagnostics = ClimaLand.default_diagnostics(model, t0, outdir), + updateat = [promote(t0:(ITime(3600 * 3)):tf...)...], + solver_kwargs = (;), ) if !isnothing(diagnostics) && !isempty(diagnostics) && + !( + first(diagnostics).output_writer isa + ClimaDiagnostics.Writers.DictWriter + ) && first(diagnostics).output_writer.output_dir != outdir @warn "Note that the kwarg outdir and outdir used in diagnostics are inconsistent; using $(first(diagnostics).output_writer.output_dir)" end - domain = ClimaLand.get_domain(model) - t0 = ITime(0, Dates.Second(1), start_date) - tf = ITime(Dates.seconds(stop_date - start_date), epoch = start_date) - Δt = ITime(Δt, epoch = start_date) - t0, tf, Δt = promote(t0, tf, Δt) # set initial conditions Y, p, cds = initialize(model) @@ -145,9 +185,10 @@ function LandSimulation( ) # Required callbacks - updateat = [promote(t0:(ITime(3600 * 3)):tf...)...] drivers = ClimaLand.get_drivers(model) updatefunc = ClimaLand.make_update_drivers(drivers) + updateat = + updateat isa AbstractVector ? updateat : [promote(t0:updateat:tf...)...] driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) required_callbacks = (driver_cb,) # TBD: can we update each step? @@ -157,16 +198,16 @@ function LandSimulation( diag_cb = ClimaDiagnostics.DiagnosticsCallback(diagnostic_handler) - # Collect all callbacks + # Collect all callbacks #TODO: ordering can be confusing as the state can be saved + # in both user_cbs and diag_cbs, and the driver update happens between them callbacks = SciMLBase.CallbackSet(user_callbacks..., required_callbacks..., diag_cb) - _integrator = SciMLBase.init( problem, timestepper; dt = Δt, callback = callbacks, - adaptive = false, + merge((; adaptive = false), solver_kwargs)..., ) return LandSimulation( model, @@ -179,6 +220,66 @@ function LandSimulation( ) end +""" + LandSimulation( + t0::AbstractFloat, + tf::AbstractFloat, + Δt::AbstractFloat, + args...; + start_date = nothing, + kwargs..., + ) + +A convenience constructor for `LandSimulation` that converts `t0`, `tf`, `Δt` into `ITime`(s), +using `start_date` as the epoch if provided. If the `kwargs` contain `updateat`, it will +convert the update times to `ITime`(s) as well. The same applies to `saveat` in `solver_kwargs`. + +""" +function LandSimulation( + t0::AbstractFloat, + tf::AbstractFloat, + Δt::AbstractFloat, + args...; + start_date = nothing, + kwargs..., +) + t0_itime, tf_itime, Δt_itime = + promote(ITime.((t0, tf, Δt); epoch = start_date)...) + kwargs = convert_kwarg_updates(t0_itime, kwargs) + LandSimulation(t0_itime, tf_itime, Δt_itime, args...; kwargs...) +end + +""" + LandSimulation( + start_date::Dates.DateTime, + stop_date::Dates.DateTime, + Δt::Union{AbstractFloat, Dates.Second}, + args...; + kwargs..., + ) + +A convenience constructor for `LandSimulation` that converts `start_date`, `stop_date`, and `Δt` +into `ITime`(s), using `start_date` as the epoch. If the `kwargs` contain `updateat`, it will +convert the update times to `ITime`(s) as well. The same applies to `saveat` in `solver_kwargs`. +""" +function LandSimulation( + start_date::Dates.DateTime, + stop_date::Dates.DateTime, + Δt::Union{AbstractFloat, Dates.Second}, + args...; + kwargs..., +) + t0 = ITime(0, Dates.Second(1), start_date) + tf = ITime( + Dates.value(convert(Dates.Second, stop_date - start_date)), + epoch = start_date, + ) + Δt = ITime(Δt isa Dates.Second ? Δt.value : Δt, epoch = start_date) + t0_itime, tf_itime, Δt_itime = promote(t0, tf, Δt) + kwargs = convert_kwarg_updates(t0_itime, kwargs) + LandSimulation(t0_itime, tf_itime, Δt_itime, args...; kwargs...) +end + """ step!(landsim::LandSimulation) @@ -236,11 +337,59 @@ model, and current date of the simulation `landsim`. function Base.show(io::IO, landsim::LandSimulation) device_type = nameof(typeof(ClimaComms.device(landsim))) model_type = nameof(typeof(landsim.model)) + t = landsim._integrator.t + final_line = if isnothing(t.epoch) + "└── Current simulation time: $(t)\n" + else + "└── Current date: $(date(t))\n" + end return print( io, "$(model_type) Simulation\n", "├── Running on: $(device_type)\n", - "└── Current date: $(date(landsim._integrator.t))\n", + final_line, ) end + +""" + convert_kwarg_updates(t0::ITime, kwargs) + +Converts the below kwargs passed to the LandSimulation constructor, which are +update vectors/intervals, into the same type as `t0`. +- `updateat`: a vector of update times, or an update interval +- `solver_kwargs[:updateat]`: a vector of update times, or an update interval +""" +function convert_kwarg_updates(t0::ITime, kwargs) + if haskey(kwargs, :solver_kwargs) + solver_kwargs = + convert_kwarg_updates(t0, kwargs[:solver_kwargs], :saveat) + kwargs = merge((; kwargs...), (; solver_kwargs,)) + end + kwargs = convert_kwarg_updates(t0, kwargs, :updateat) + return kwargs +end + +function convert_kwarg_updates(t0::ITime, kwargs, key::Symbol) + haskey(kwargs, key) || return kwargs # if the key is not present, return the original kwargs + updates = kwargs[key] + converted_updates = convert_updates.(t0, updates) + if converted_updates isa Tuple + converted_updates = collect(converted_updates) + end + return merge((; kwargs...), (; key => converted_updates)) +end + +convert_updates(t0::ITime, update_time::ITime) = promote(t0, update_time)[2] +convert_updates(t0::ITime, update_time::Dates.DateTime) = promote( + t0, + ITime( + Dates.value(convert(Dates.Second, update_time - t0.epoch)); + epoch = t0.epoch, + ), +)[2] +convert_updates(t0::ITime, update_time::Dates.Period) = + promote(t0, ITime(Dates.value(convert(Dates.Second, update_time));))[2] +convert_updates(t0::ITime, update_time::AbstractFloat) = + promote(t0, ITime(update_time, epoch = t0.epoch))[2] + end#module diff --git a/src/standalone/Soil/energy_hydrology.jl b/src/standalone/Soil/energy_hydrology.jl index 4d5cf07695..47f7e3413a 100644 --- a/src/standalone/Soil/energy_hydrology.jl +++ b/src/standalone/Soil/energy_hydrology.jl @@ -534,9 +534,11 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT} @. p.soil.full_bidiag_matrix_scratch += MatrixFields.LowerDiagonalMatrixRow(p.soil.topBC_scratch) end - + # dtγ can be an ITime or a float @. ∂ϑres∂ϑ = - -dtγ * (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) + FT(-1) * + float(dtγ) * + (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) # Now create the flux term for ∂ρe∂ϑ using bidiag_matrix_scratch # This overwrites full_bidiag_matrix_scratch @@ -550,7 +552,9 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT} ), ) ⋅ p.soil.bidiag_matrix_scratch @. ∂ρeres∂ϑ = - -dtγ * (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) + FT(-1) * + float(dtγ) * + (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) # Now overwrite bidiag_matrix_scratch and full_bidiag scratch for the ρe ρe bidiagonal @. p.soil.bidiag_matrix_scratch = @@ -566,7 +570,9 @@ function ClimaLand.make_compute_jacobian(model::EnergyHydrology{FT}) where {FT} MatrixFields.DiagonalMatrixRow(interpc2f_op(-p.soil.κ)) ⋅ p.soil.bidiag_matrix_scratch @. ∂ρeres∂ρe = - -dtγ * (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) + FT(-1) * + float(dtγ) * + (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) end return compute_jacobian! end diff --git a/src/standalone/Soil/rre.jl b/src/standalone/Soil/rre.jl index f33a1da6ed..f99a9e2e9f 100644 --- a/src/standalone/Soil/rre.jl +++ b/src/standalone/Soil/rre.jl @@ -438,7 +438,9 @@ function ClimaLand.make_compute_jacobian(model::RichardsModel{FT}) where {FT} end # Compute divergence matrix @. ∂ϑres∂ϑ = - -dtγ * (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) + FT(-1) * + float(dtγ) * + (divf2c_matrix() ⋅ p.soil.full_bidiag_matrix_scratch) - (I,) end diff --git a/src/standalone/Vegetation/canopy_energy.jl b/src/standalone/Vegetation/canopy_energy.jl index ca83005712..fdf04ba1a0 100644 --- a/src/standalone/Vegetation/canopy_energy.jl +++ b/src/standalone/Vegetation/canopy_energy.jl @@ -190,7 +190,7 @@ function ClimaLand.make_compute_jacobian( Y.canopy.energy.T - _T_freeze, )# use atmos air pressure as approximation for surface air pressure @. ∂Tres∂T = - dtγ * MatrixFields.DiagonalMatrixRow( + float(dtγ) * MatrixFields.DiagonalMatrixRow( (∂LW_n∂Tc - ∂SHF∂Tc - ∂LHF∂qc * ∂qc∂Tc) / (ac_canopy * max(area_index.leaf + area_index.stem, eps(FT))), ) - (I,) @@ -241,7 +241,7 @@ end A function which updates `surface_field` in place with the value of the big leaf model's energy. -Note that this assumes that there is at most a single canopy and stem +Note that this assumes that there is at most a single canopy and stem component, and that the area index for them refers to the integrated area index (in height) - not the value per layer. """ diff --git a/test/shared_utilities/utilities.jl b/test/shared_utilities/utilities.jl index 239184d1e4..1e55b7538d 100644 --- a/test/shared_utilities/utilities.jl +++ b/test/shared_utilities/utilities.jl @@ -4,7 +4,8 @@ ClimaComms.@import_required_backends using ClimaCore: Spaces, Geometry, Fields using ClimaLand using ClimaLand: Domains, condition, SavingAffect, saving_initialize - +using Dates +import ClimaUtilities.TimeManager: ITime, date ## Callback tests mutable struct Integrator{FT} t::Any @@ -104,6 +105,66 @@ end @test cb.affect!.saved_values.saveval[1] == integrator.p end +@testset "NonInterpSavingCallback" begin + start_date = ITime(0, Dates.Second(1), DateTime(2020)) + stop_date = start_date + ITime(60.0 * 60) + dt = ITime(60.0) + save_interval = ITime(60.0 * 20) + + saveat_itime = collect(start_date:save_interval:stop_date) + sv_itime = (; + t = Array{Union{Nothing, eltype(saveat_itime)}}( + nothing, + length(saveat_itime), + ), + saveval = Array{Any}(undef, length(saveat_itime)), + ) + saving_cb_itime = ClimaLand.NonInterpSavingCallback(sv_itime, saveat_itime) + + saveat_ft = collect(0:float(save_interval):float(stop_date - start_date)) + sv_ft = (; + t = Array{Union{Nothing, eltype(saveat_ft)}}( + nothing, + length(saveat_ft), + ), + saveval = Array{Any}(undef, length(saveat_ft)), + ) + saving_cb_ft = ClimaLand.NonInterpSavingCallback(sv_ft, saveat_ft) + + saveat_date = collect( + date(start_date):Dates.Second(float(save_interval)):date(stop_date), + ) + sv_date = (; + t = Array{Union{Nothing, eltype(saveat_date)}}( + nothing, + length(saveat_date), + ), + saveval = Array{Any}(undef, length(saveat_date)), + ) + saving_cb_date = ClimaLand.NonInterpSavingCallback(sv_date, saveat_date) + + for t in start_date:dt:stop_date + saving_cb_itime.condition(nothing, t, nothing) && + saving_cb_itime.affect!(Integrator(t, [t])) + saving_cb_ft.condition(nothing, t, nothing) && + saving_cb_ft.affect!(Integrator(t, [float(t)])) + saving_cb_date.condition(nothing, t, nothing) && + saving_cb_date.affect!(Integrator(t, [date(t)])) + end + for (i, t) in enumerate(saveat_itime) + @test sv_itime.t[i] == t + @test sv_itime.saveval[i] == [t] + end + for (i, t) in enumerate(saveat_ft) + @test sv_ft.t[i] == t + @test sv_ft.saveval[i] == [t] + end + for (i, t) in enumerate(saveat_date) + @test sv_date.t[i] == t + @test sv_date.saveval[i] == [t] + end +end + ## DSS 0D/1D tests @testset "dss! - 0D/1D spaces, FT = $FT" begin