diff --git a/CHANGELOG.md b/CHANGELOG.md index 1975b1a..0b330ae 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,7 @@ +# v1.4 + +- Moved to MTKv11 (ProcessBasedModelling.jl v1.9) + # v1.3 ## `CloudToppedMixedLayerModel` diff --git a/Project.toml b/Project.toml index f5137ec..90964b8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,12 @@ name = "ConceptualClimateModels" uuid = "3e640dcb-3446-4ed7-aadb-0366b77312ec" authors = ["Datseris "] -version = "1.3.0" +version = "1.4.0" [deps] DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4" NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" +OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" ProcessBasedModelling = "ca969041-2cf3-4b10-bc21-86f4417093eb" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" @@ -13,7 +14,8 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" [compat] DynamicalSystemsBase = "3.5" NaNMath = "1" -ProcessBasedModelling = "1.8" +OrdinaryDiffEqNonlinearSolve = "2" +ProcessBasedModelling = "1.9" Reexport = "1.2" Roots = "2" julia = "1.9" \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 5d43c27..2813402 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,6 +34,6 @@ bib = CitationBibliography( build_docs_with_style(pages, ConceptualClimateModels, GlobalMeanEBM, CloudToppedMixedLayerModel, ProcessBasedModelling; authors = "George Datseris ", - bib, warnonly = [:doctest, :missing_docs, :cross_references], + bib, warnonly = [:doctest, :missing_docs, :cross_references, :linkcheck], repo = Remotes.GitHub("JuliaDynamics", "ConceptualClimateModels.jl") ) diff --git a/src/dynamicalsystems.jl b/src/dynamicalsystems.jl index c7ccf2b..b8d8790 100644 --- a/src/dynamicalsystems.jl +++ b/src/dynamicalsystems.jl @@ -7,33 +7,42 @@ export all_equations export named_current_parameters export try_set_parameter! +# This is needed because otherwise we get the error: +# This ODE requires a DAE initialization and thus a nonlinear solve but no nonlinear solve +# has been loaded. To solve this problem, do `using OrdinaryDiffEqNonlinearSolve` or pass +# a custom `nlsolve` choice into the `initializealg`. +using OrdinaryDiffEqNonlinearSolve + DEFAULT_DIFFEQ = DynamicalSystemsBase.DEFAULT_DIFFEQ """ processes_to_coupledodes(processes [, default]; kw...) -Convert a given `Vector` of processes to a `DynamicalSystem`, in particular `CoupledODEs`. +Convert a given `Vector` of processes to `CoupledODEs`. All processes represent symbolic equations managed by ModelingToolkit.jl. `default` is a vector for default processes that "process-less" variables introduced in `processes` will obtain. -Use [`processes_to_mtkmodel`](@ref) to obtain the MTK model before it is structurally -simplified and converted to a `DynamicalSystem`. +Use [`processes_to_mtkmodel`](@ref) to obtain the MTK model before it is compiled +and converted to `CoupledODEs`. See also [`processes_to_mtkmodel`](@ref) for more details on what `processes` is, or see the online [Tutorial](@ref). ## Keyword arguments -- `diffeq`: options passed to DifferentialEquations.jl ODE solving +- `diffeq`: options passed to OrdinaryDiffEq.jl ODE solving when constructing the `CoupledODEs`. - `inplace`: whether the dynamical system will be in place or not. Defaults to `true` if the system dimension is ≤ 5. - `split = false`: whether to split parameters as per ModelingToolkit.jl. Note the default is not ModelingToolkit's default, i.e., no splitting occurs. - This accelerates parameter access, assuming all parameters are of the same type. + This accelerates parameter access by assuming all parameters are of the same type. +- `probkw = (warn_initialize_determined = false, )`: keyword arguments expanded into + the creation of `ODEProblem`. - `kw...`: all other keywords are propagated to `processes_to_mtkmodel`. """ function processes_to_coupledodes(proc, default = []; - diffeq = DEFAULT_DIFFEQ, inplace::Bool = false, split::Bool = false, kwargs... + diffeq = DEFAULT_DIFFEQ, inplace = nothing, split::Bool = false, + probkw = (warn_initialize_determined = false, ), kwargs... ) sys = processes_to_mtkmodel(proc, default; kwargs...) ssys = mtkcompile(sys; split) @@ -45,9 +54,9 @@ function processes_to_coupledodes(proc, default = []; # The usage of `nothing` for the initial state assumes all state variables # and all parameters have been defined with a default value. This also means if inplace - prob = ODEProblem(ssys, nothing, (0.0, Inf)) + prob = ODEProblem(ssys, nothing, (0.0, Inf); probkw...) else - prob = ODEProblem{false}(ssys, nothing, (0.0, Inf); u0_constructor = x->SVector(x...)) + prob = ODEProblem{false}(ssys, nothing, (0.0, Inf); u0_constructor = x->SVector(x...), probkw...) end ds = CoupledODEs(prob, diffeq) return ds diff --git a/test/Project.toml b/test/Project.toml index fb09431..7268173 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,5 @@ [deps] -DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" +Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/test/runtests.jl b/test/runtests.jl index 6343f44..651bec0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,88 +3,85 @@ using ConceptualClimateModels.GlobalMeanEBM using Test @testset "all processes" begin -function test_symbolic_var(mtk, var) - if has_symbolic_var(mtk, var) - @test true - else - error("var $var doesn't exist") + function test_symbolic_var(mtk, var) + if has_symbolic_var(mtk, var) + @test true + else + error("var $var doesn't exist") + end end -end - -@variables begin - ε1(t) = 0.1 - ε2(t) = 0.1 - ε3(t) = 0.1 - ε4(t) = 0.1 -end -p1 = [ - BasicRadiationBalance(), - EmissivityStefanBoltzmanOLR(), - IceAlbedoFeedback(; min = 0.3, max = 0.7), - α ~ α_ice, - ParameterProcess(ε), # emissivity is a parameter - f ~ 0, # no external forcing - # absorbed solar radiation has a default process - EmissivityFeedbackTanh(ε = ε1), - SoedergrenClearSkyEmissivity(ε = ε2), - EmissivitySellers1969(ε = ε3), - EmissivityLogConcentration(ε = ε4), -] - -ds = processes_to_coupledodes(p1, GlobalMeanEBM) -mtk = referrenced_sciml_model(ds) -eqs = all_equations(mtk) - -# This also tests a bunch of default processes -for var in (T, ε, ε1, ε2, ε3, ε4, q, :ε1_T_sigmoid_ref, mtk.ε_0) - test_symbolic_var(eqs, var) -end + @variables begin + ε1(t) = 0.1 + ε2(t) = 0.1 + ε3(t) = 0.1 + ε4(t) = 0.1 + end + p1 = [ + BasicRadiationBalance(), + EmissivityStefanBoltzmanOLR(), + IceAlbedoFeedback(; min = 0.3, max = 0.7), + α ~ α_ice, + ParameterProcess(ε), # emissivity is a parameter + f ~ 0, # no external forcing + # absorbed solar radiation has a default process + EmissivityFeedbackTanh(ε = ε1), + SoedergrenClearSkyEmissivity(ε = ε2), + EmissivitySellers1969(ε = ε3), + EmissivityLogConcentration(ε = ε4), + ] -@variables begin - OLR1(t) = 10.0 - OLR2(t) = 10.0 - OLR3(t) = 10.0 - a1(t) = 0.1 - a2(t) = 0.1 - a3(t) = 0.1 - a4(t) = 0.1 - ΔS1(t) = 0.0 - ΔT1(t) = 0.0 -end + ds = processes_to_coupledodes(p1, GlobalMeanEBM) + mtk = referrenced_sciml_model(ds) + eqs = all_equations(mtk) -p2 = [ - LinearOLR(), - LinearClearSkyOLR(; OLR = OLR1), - BudykoOLR(; OLR = OLR2), - CloudAlbedoExponential(), - CloudAlbedoLinear(; α_cloud = a1), - DirectAlbedoAddition(α = a4), - CoAlbedoProduct(α = a2), - SeparatedClearAllSkyAlbedo(α = a3), - C ~ 0.6, - ΔTLinearRelaxation(), - ΔTStommelModel(ΔT = ΔT1, ΔS = ΔS1), -] - -ds = processes_to_coupledodes(p2, GlobalMeanEBM) -mtk = referrenced_sciml_model(ds) -for var in (T, C, OLR1, OLR2, :α_bg, a1, a2, a3, a4, ΔS1, ΔT) - if has_symbolic_var(mtk, var) - @test true - else - error("var $var doesn't exist") + # This also tests a bunch of default processes + for var in (T, ε, ε1, ε2, ε3, ε4, q, :ε1_T_sigmoid_ref, mtk.ε_0) + test_symbolic_var(eqs, var) end -end + @variables begin + OLR1(t) = 10.0 + OLR2(t) = 10.0 + OLR3(t) = 10.0 + a1(t) = 0.1 + a2(t) = 0.1 + a3(t) = 0.1 + a4(t) = 0.1 + ΔS1(t) = 0.0 + ΔT1(t) = 0.0 + end + + p2 = [ + LinearOLR(), + LinearClearSkyOLR(; OLR = OLR1), + BudykoOLR(; OLR = OLR2), + CloudAlbedoExponential(), + CloudAlbedoLinear(; α_cloud = a1), + DirectAlbedoAddition(α = a4), + CoAlbedoProduct(α = a2), + SeparatedClearAllSkyAlbedo(α = a3), + C ~ 0.6, + ΔTLinearRelaxation(), + ΔTStommelModel(ΔT = ΔT1, ΔS = ΔS1), + ] + ds = processes_to_coupledodes(p2, GlobalMeanEBM) + mtk = referrenced_sciml_model(ds) + for var in (T, C, OLR1, OLR2, :α_bg, a1, a2, a3, a4, ΔS1, ΔT) + if has_symbolic_var(mtk, var) + @test true + else + error("var $var doesn't exist") + end + end end @testset "dynamical systems integration" begin - using DynamicalSystems + using Attractors budyko_processes = [ BasicRadiationBalance(), @@ -100,7 +97,6 @@ end grid = plausible_grid(budyko) mapper = AttractorsViaRecurrences(budyko, grid) - rfam = RecurrencesFindAndMatch(mapper) sampler = plausible_ic_sampler(budyko) set_parameter!(budyko, :ε_0, 0.5) @@ -109,5 +105,4 @@ end attractors = extract_attractors(mapper) @test length(attractors) == 2 - end \ No newline at end of file