Merge remote-tracking branch 'origin/master' into fix_ode_jumps

Author: isaacsas

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Yingbo Ma <[email protected]>", "Chris Rackauckas <[email protected]> and contributors"]
-version = "9.68.0"
+version = "9.68.1"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

docs/Project.toml

@@ -29,7 +29,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 [compat]
 BenchmarkTools = "1.3"
 BifurcationKit = "0.4"
-DataInterpolations = "6.5"
+DataInterpolations = "6.5, 8"
 Distributions = "0.25"
 Documenter = "1"
 DynamicQuantities = "^0.11.2, 0.12, 1"

docs/src/comparison.md

@@ -23,9 +23,9 @@
   - Modelica is an object-oriented single dispatch language. ModelingToolkit.jl,
     built on Julia, uses multiple dispatch extensively to simplify code.
   - Many Modelica compilers supply a GUI. ModelingToolkit.jl does not.
-  - Modelica can be used to simulate ODE and DAE systems. ModelingToolkit.jl
-    has a much more expansive set of system types, including nonlinear systems,
-    SDEs, PDEs, and more.
+  - Modelica is designed for simulating ODE and DAE systems (which can include nonlinear dynamics).
+    In contrast, ModelingToolkit.jl supports a much broader range of system types, including SDEs,
+    PDEs, time-independent nonlinear systems (e.g. various forms of optimization problems) and more.
 
 ## Comparison Against Simulink
 

docs/src/tutorials/change_independent_variable.md

@@ -4,9 +4,10 @@ Ordinary differential equations describe the rate of change of some dependent va
 For the modeler it is often most natural to write down the equations with a particular independent variable, say time $t$.
 However, in many cases there are good reasons for changing the independent variable:
 
-    1. One may want $y(x)$ as a function of $x$ instead of $(x(t), y(t))$ as a function of $t$
-    2. Some differential equations vary more nicely (e.g. less stiff) with respect to one independent variable than another.
-    3. It can reduce the number of equations that must be solved (e.g. $y(x)$ is one equation, while $(x(t), y(t))$ are two).
+ 1. One may want $y(x)$ as a function of $x$ instead of $(x(t), y(t))$ as a function of $t$
+
+ 2. Some differential equations vary more nicely (e.g. less stiff) with respect to one independent variable than another.
+ 3. It can reduce the number of equations that must be solved (e.g. $y(x)$ is one equation, while $(x(t), y(t))$ are two).
 
 To manually change the independent variable of an ODE, one must rewrite all equations in terms of a new variable and transform differentials with the chain rule.
 This is mechanical and error-prone.
@@ -34,9 +35,9 @@ It expresses the position $(x(t), y(t))$ as a function of time $t$.
 But suppose we want to determine whether the projectile hits a target 10 meters away.
 There are at least three ways of answering this:
 
-    - Solve the ODE for $(x(t), y(t))$ and use a callback to terminate when $x$ reaches 10 meters, and evaluate $y$ at the final time.
-    - Solve the ODE for $(x(t), y(t))$ and use root finding to find the time when $x$ reaches 10 meters, and evaluate $y$ at that time.
-    - Solve the ODE for $y(x)$ and evaluate it at 10 meters.
+  - Solve the ODE for $(x(t), y(t))$ and use a callback to terminate when $x$ reaches 10 meters, and evaluate $y$ at the final time.
+  - Solve the ODE for $(x(t), y(t))$ and use root finding to find the time when $x$ reaches 10 meters, and evaluate $y$ at that time.
+  - Solve the ODE for $y(x)$ and evaluate it at 10 meters.
 
 We will demonstrate the last method by changing the independent variable from $t$ to $x$.
 This transformation is well-defined for any non-zero horizontal velocity $v$, so $x$ and $t$ are one-to-one.
@@ -55,7 +56,7 @@ M2s # display this # hide
 
 The derivatives are now with respect to the new independent variable $x$, which can be accessed with `M2.x`.
 
-!!! warn
+!!! info
 
     At this point `x`, `M1.x`, `M1s.x`, `M2.x`, `M2s.x` are *three* different variables.
     Meanwhile `y`, `M1.y`, `M1s.y`, `M2.y` and `M2s.y` are *four* different variables.
@@ -66,7 +67,7 @@ It is straightforward to evolve the ODE for 10 meters and plot the resulting tra
 
 ```@example changeivar
 using OrdinaryDiffEq, Plots
-prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters with x-velocity 8 m/s
+prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters
 sol = solve(prob, Tsit5())
 plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)!
 ```

src/systems/jumps/jumpsystem.jl

@@ -4,25 +4,14 @@ const JumpType = Union{VariableRateJump, ConstantRateJump, MassActionJump}
 # call reset_aggregated_jumps!(integrator).
 # assumes iip
 function _reset_aggregator!(expr, integrator)
-    if expr isa Symbol
-        error("Error, encountered a symbol. This should not happen.")
+    @assert Meta.isexpr(expr, :function)
+    body = expr.args[end]
+    body = quote
+        $body
+        $reset_aggregated_jumps!($integrator)
     end
-    if expr isa LineNumberNode
-        return
-    end
-
-    if (expr.head == :function)
-        _reset_aggregator!(expr.args[end], integrator)
-    else
-        if expr.args[end] == :nothing
-            expr.args[end] = :(reset_aggregated_jumps!($integrator))
-            push!(expr.args, :nothing)
-        else
-            _reset_aggregator!(expr.args[end], integrator)
-        end
-    end
-
-    nothing
+    expr.args[end] = body
+    return nothing
 end
 
 """

src/systems/parameter_buffer.jl

@@ -444,6 +444,9 @@ end
 
 function validate_parameter_type(ic::IndexCache, stype, sz, sym, index, val)
     (; portion) = index
+    if stype <: FnType
+        stype = fntype_to_function_type(stype)
+    end
     # Nonnumeric parameters have to match the type
     if portion === NONNUMERIC_PORTION
         val isa stype && return nothing

test/initial_values.jl

@@ -1,6 +1,7 @@
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D, get_u0
 using OrdinaryDiffEq
+using DataInterpolations
 using SymbolicIndexingInterface: getu
 
 @variables x(t)[1:3]=[1.0, 2.0, 3.0] y(t) z(t)[1:2]
@@ -219,3 +220,19 @@ end
     @test_throws ["a(t)", "c(t)"] ODEProblem(
         sys, [e => 2, a => b, b => a + 1, c => d, d => c + 1], (0, 1))
 end
+
+@testset "Issue#3490: `remake` works with callable parameters" begin
+    ts = collect(0.0:0.1:10.0)
+    spline = LinearInterpolation(ts .^ 2, ts)
+    Tspline = typeof(spline)
+    @variables x(t)
+    @parameters (interp::Tspline)(..)
+
+    @mtkbuild sys = ODESystem(D(x) ~ interp(t), t)
+
+    prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), [interp => spline])
+    spline2 = LinearInterpolation(ts .^ 2, ts .^ 2)
+    p_new = [interp => spline2]
+    prob2 = remake(prob; p = p_new)
+    @test prob2.ps[interp] == spline2
+end