modified: .gitignore /.vscode

modified:   Project.toml Version, extras
modified:   README.md    First version
modified:   src/MechGluecode.jl  export norm
new file:   test/diffeq_1.jl    First version
new file:   test/diffeq_plots_1.jl First version
new file:   test/diffeq_plots_2.jl First version
new file:   test/plots_1.jl First version
modified:   test/runtests.jl Include all, comment snoopcompile

.gitignore

@@ -1 +1,2 @@
 /Manifest.toml
+/.vscode

Project.toml

@@ -1,7 +1,7 @@
 name = "MechGluecode"
 uuid = "3017d99d-ab52-4519-99a2-fa9ddc4637fe"
 authors = ["hustf <hustf@users.noreply.github.com> and contributors"]
-version = "0.1.4"
+version = "0.1.5"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -24,6 +24,8 @@ julia = "1"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+SnoopCompile = "aa65fe97-06da-5843-b5b1-d5d13cad87d2"
+SnoopCompileCore = "e2b509da-e806-4183-be48-004708413034"
 
 [targets]
 test = ["Test"]

README.md

@@ -1 +1,13 @@
 # MechGluecode
+
+Glue code for the types defined in [Unitfu.jl](https://github.com/hustf/Unitfu.jl). The required packages are registered in [M8](https://github.com/hustf/M8).
+
+The extension glue code is loaded if you also load:
+    - DifferentialEquations.jl => load MechGlueDiffEqBase.jl
+    - Plots.jl => load MechGluePlots.jl
+
+Not completed yet:
+    - Interpolations.jl => load MechGlueDiffEqBase.jl
+    - Plots.jl => load MechGlueInterpolations.jl
+    - ModelingToolkit => MechGlueModelingToolkit
+    - [RecursiveArrayTools](https://github.com/SciML/RecursiveArrayTools.jl): This is included in MechGlueDiffEqBase, no standalone glue code needed (?)

src/MechGluecode.jl

@@ -1,37 +1,37 @@
 module MechGluecode
 using Requires
-export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu
+export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu, norm
 
 function __init__()
     @require MechanicalUnits = "e6be9192-89dc-11e9-36e6-5dbcb28f419e" begin
         @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
             @info "Plots => using MechGluePlots"
             @eval using MechGluePlots
         end
-        @require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" begin
-            import MechanicalUnits: Unitfu, @import_expand
-            import Unitfu: AbstractQuantity
-            import DifferentialEquations
-            import DifferentialEquations: DiffEqBase
-            import DiffEqBase: value, ODE_DEFAULT_NORM, UNITLESS_ABS2
-            @info "DifferentialEquations => using MechGlueDiffEqBase"
-            @eval using MechGlueDiffEqBase
-        end
-
         @require Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" begin
-            @info "Interpolations => using MechGlueInterpolations"
-            @eval using MechGlueInterpolations
+            #@info "Interpolations => using MechGlueInterpolations"
+            #@eval using MechGlueInterpolations
         end
-        
+
         @require ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" begin
-            @info "ModelingToolkit => using MechGlueModelingToolkit"
-            @eval using MechGlueModelingToolkit
+            #@info "ModelingToolkit => using MechGlueModelingToolkit"
+            #@eval using MechGlueModelingToolkit
         end
 
         @require RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" begin
-            @eval using MechGlueDiffEqBase
             @info "RecursiveArrayTools => using MechGlueRecursiveArrayTools"
-            @eval using MechGlueRecursiveArrayTools
+           # @eval using MechGlueDiffEqBase
+           #  @info "RecursiveArrayTools => using MechGlueRecursiveArrayTools"
+           #  @eval using MechGlueRecursiveArrayTools
+        end
+        @require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" begin
+            @info "DifferentialEquations => using MechGlueDiffEqBase"
+          #  @eval import MechanicalUnits: Unitfu, @import_expand
+          #  @eval import Unitfu: AbstractQuantity
+          #  @eval import DifferentialEquations
+          #  @eval import DifferentialEquations: DiffEqBase
+          #  @eval import SciMLBase
+            @eval using MechGlueDiffEqBase
         end
     end
     @info "MechGluecode init"

test/diffeq_1.jl

@@ -0,0 +1,28 @@
+using Test
+using MechGluecode
+using MechanicalUnits
+using MechanicalUnits: Unitfu.promote_to_derived
+@import_expand(m, s, N) # Takes long if defined late)
+using DifferentialEquations
+
+
+@testset "Unitful time with unitless state" begin
+    u0 = 30.0s
+    tspan = (0.0, 10.0)s
+    prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
+    prob2 = ODEProblem((u,t,p)    -> (-0.2∙s⁻¹ * u[1]), u0, tspan)
+    prob3 = ODEProblem((u,t,p)    -> [-0.2∙s⁻¹* u[1]], [u0],tspan)
+    for prob in [prob1, prob2, prob3]
+        @test solve(prob, Tsit5()).retcode === :Success
+    end
+end
+@testset "ODE quantity" begin
+    f = (y, p, t) -> 0.5y / 3.0s  # The derivative, yeah
+    u0 = 1.5N
+    tspan = (0.0s, 1.0s)
+    prob = ODEProblem(f, u0, tspan)
+    sol = solve(prob)
+    @test sol.t isa Vector{<:Time}
+    @test sol.u isa Vector{<:Force}
+end
+

test/diffeq_plots_1.jl

@@ -0,0 +1,120 @@
+using Test
+using MechGluecode
+using MechanicalUnits
+using MechanicalUnits: Unitfu.promote_to_derived
+@import_expand(m, s, N) # Takes long if defined late)
+using DifferentialEquations
+
+
+@testset "Unitful time with unitful function" begin
+    u0 = 30.0N
+    tspan = (0.0, 10.0)s
+    prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
+    prob2 = ODEProblem((u,t,p)    ->         (-0.2∙s⁻¹ * u[1]), u0, tspan)
+    prob3 = ODEProblem((u,t,p)    ->         [-0.2∙s⁻¹* u[1]], [u0],tspan)
+    so1 =  solve(prob1, Euler(), dt = 1s)
+    so2 =  solve(prob2, Midpoint())
+    so3 =  solve(prob3, Tsit5())
+    plot(so1, denseplot = false, marker = true)
+    plot(so2, denseplot = false, marker = true)
+    plot(so3, denseplot = false, marker = true)
+    sf1(t) = so1(t)[1]
+    sf2(t) = so2(t)[1]
+    sf3(t) = so3(t)[1]
+    plot([sf1, sf2, sf3], xlims =(0,10)s, ylims =(0,40)N)
+    plot([t-> sf1(t) / sf3(t), t-> sf2(t) / sf3(t)], xlims =(0,10)s, ylims = (0.5,1.1), minorgrid=true, grid = (:y, :olivedrab, :dot, 1, 0.9))
+
+end
+@testset "Unitful time with unitless function" begin
+    u0 = 30.0
+    tspan = (0.0, 10.0)s
+    prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
+    prob2 = ODEProblem((u,t,p)    ->         (-0.2∙s⁻¹ * u[1]), u0, tspan)
+    prob3 = ODEProblem((u,t,p)    ->         [-0.2∙s⁻¹* u[1]], [u0],tspan)
+    so1 =  solve(prob1, Euler(), dt = 1s)
+    so2 =  solve(prob2, Midpoint())
+    so3 =  solve(prob3, Tsit5())
+    # don't work plot(so1, denseplot = false, marker = true)
+    # don't work plot(so2, denseplot = false, marker = true)
+    # don't work plot(so3, denseplot = false, marker = true)
+    sf1(t) = so1(t)[1]
+    sf2(t) = so2(t)[1]
+    sf3(t) = so3(t)[1]
+    plot([sf1, sf2, sf3], xlims =(0,10)s, ylims =(0,40))
+    plot([t-> sf1(t) / sf3(t), t-> sf2(t) / sf3(t)], xlims =(0,10)s, ylims = (0.5,1.1), minorgrid=true, grid = (:y, :olivedrab, :dot, 1, 0.9))
+end
+
+@testset "Plot ODE quantity" begin
+    # Don't give us kN from from solution in this case, just defaults.
+    promote_to_derived()
+    f1 = (y, p, t) -> -1.5y / 0.3s
+    u0 = 1.5N
+    tspan = (0.0s, 1.0s)
+    prob = ODEProblem(f1, u0, tspan)
+    sol = solve(prob)
+    plot(sol)
+end
+
+@testset "Plot ODE system same quantities" begin
+    # This is not a good example, as the units are incorrect and have no
+    # known physical interpretation
+    σ = 10/s
+    ρ = 28m
+    β = (8 / 3)m
+    function lorenz!(du, u, p, t)
+        du[1] = (u[2] - u[1])∙σ
+        du[2] = (u[1] * (ρ - u[3])) / (m∙s) - u[2]/s
+        du[3] =  u[1] * u[2] / (m∙s)    - β ∙ u[3]/(m∙s)
+    end
+    u0 = [1.0, 0.0, 0.0]m
+    tspan = (0.0,100.0)s
+    prob = ODEProblem(lorenz!,u0,tspan)
+    sol = solve(prob)
+    # Plot all variables vs time
+    plot(sol)
+    # Plot in phase space
+    plot(sol, vars=(1,2,3))
+end
+#=
+@testset "Plot ODE system single pendulum" begin
+    # Second order non-linear, undamped pendulum:
+    # mLθ'' + mg∙sin(θ)
+    # where m is mass, unit kg
+    #       L is length, unit m
+    #       θ is angle to vertical, strictly no unit
+    #       g is gravity's acceleration
+    # Acceleration
+    # θ'' + k * sin(θ) = 0
+    #    where k = g / L , unit(k) = s⁻²
+    # Reduced to first order non-linear system by defining:
+    # θ = y[1]
+    # θ' = y[2]
+    # It follows that
+    # θ'' = y'[2]
+    # and
+    # y'[1] = y[2]
+    # y'[2] = -k * sin(y[1])
+    # To solve numerically, we define the derivative.
+    # For efficiency, it's evaluated in place, i.e. dy changes at every call.
+    function f!(dy, y, k, t)
+        dy[1] = y[2]
+        dy[2] = -k ∙ sin(y[1])
+    end
+    # Fixed parameter
+    k = upreferred(g / 1.0m)
+    # Initial conditions
+    θ0 = 0.1
+    θ0_vel = 0.0s⁻¹
+    θ0_acc = -k *sin(θ0)
+
+    y0 = [θ0, θ0_vel]
+    y0 = ArrayPartition([θ0_vel, θ0_acc])
+    tspan = (0.0, 100.0)s
+    prob = ODEProblem(f!, y0, tspan, k)
+    sol = solve(prob)
+    # Plot all variables vs time
+    plot(sol)
+    # Plot in phase space
+    plot(solve(prob), vars=(1,2))
+end
+=#

test/diffeq_plots_2.jl

@@ -0,0 +1,84 @@
+using Test
+using MechGluecode
+using MechanicalUnits
+using MechanicalUnits: Unitfu.promote_to_derived
+@import_expand(~m, s, N)
+using DifferentialEquations
+import MechGluecode: norm # should work without
+
+@testset "Initial checks ArrayPartition" begin
+    r0 = [1131.340, -2282.343, 6672.423]∙km
+    v0 = [-5.64305, 4.30333, 2.42879]∙km/s
+    Δt = 86400.0*365∙s
+    μ = 398600.4418∙km³/s²
+    rv0 = ArrayPartition(r0,v0)
+
+    function fo(dy, y, μ, t)
+        r = norm(y.x[1])
+        dy.x[1] .= y.x[2]
+        dy.x[2] .= -μ .* y.x[1] / r^3
+    end
+    prob = ODEProblem(fo, rv0, (0.0s, 0.035Δt), μ)
+    @time sol = solve(prob, Tsit5(), dt = 0.1s);
+    # Plot all variables vs time
+    plot(sol)
+    plot(sol, vars = (1))
+    plot(sol, vars = (2))
+    plot(sol, vars = (3))
+    plot(sol, vars = (4))
+    plot(sol, vars = (5))
+    plot(sol, vars = (6))
+
+    # Plot in phase space
+    plot(sol, vars=(1,4))
+    plot(sol, vars=(2,5))
+    plot(sol, vars=(3,6))
+    plot(sol, vars=(2,5))
+    # Three distances vs time (0)
+    plot(sol, vars=[(0,1), (0,2),(0,3)])
+    # Three velocities vs time (0)
+    plot(sol, vars=[(0,4), (0,5),(0,6)])
+end
+#=
+@testset "Plot ODE system single pendulum" begin
+    # Second order non-linear, undamped pendulum:
+    # mLθ'' + mg∙sin(θ)
+    # where m is mass, unit kg
+    #       L is length, unit m
+    #       θ is angle to vertical, strictly no unit
+    #       g is gravity's acceleration
+    # Acceleration
+    # θ'' + k * sin(θ) = 0
+    #    where k = g / L , unit(k) = s⁻²
+    # Reduced to first order non-linear system by defining:
+    # θ = y[1]
+    # θ' = y[2]
+    # It follows that
+    # θ'' = y'[2]
+    # and
+    # y'[1] = y[2]
+    # y'[2] = -k * sin(y[1])
+    # To solve numerically, we define the derivative.
+    # For efficiency, it's evaluated in place, i.e. dy changes at every call.
+    function f!(dy, y, k, t)
+        dy[1] = y[2]
+        dy[2] = -k ∙ sin(y[1])
+    end
+    # Fixed parameter
+    k = upreferred(g / 1.0m)
+    # Initial conditions
+    θ0 = 0.1
+    θ0_vel = 0.0s⁻¹
+    θ0_acc = -k *sin(θ0)
+
+    y0 = [θ0, θ0_vel]
+    y0 = ArrayPartition([θ0_vel, θ0_acc])
+    tspan = (0.0, 100.0)s
+    prob = ODEProblem(f!, y0, tspan, k)
+    sol = solve(prob)
+    # Plot all variables vs time
+    plot(sol)
+    # Plot in phase space
+    plot(solve(prob), vars=(1,2))
+end
+=#

test/plots_1.jl

@@ -0,0 +1,30 @@
+using Test
+import MechanicalUnits: @import_expand
+@import_expand(m, s, N) # Takes long if defined late)
+using MechGluecode
+import Plots
+import MechGluecode
+
+@testset "Plots mixed quantity" begin
+    s1x = [1.0,2,3]
+    s2x = [1.0,2,3.5]
+    s1y = [1.0,2,4]
+    s2y = [1.0,2,2]
+    mxqm = hcat(s1x∙m, s2x∙m²)
+    myqm = hcat(s1y∙s, s2y∙s²)
+    plot(mxqm, myqm, seriestype = [:path :scatter], label = ["path" "scatter"], xguide = "x", yguide = "y")
+end
+
+@testset "Plots functions with mixed units" begin
+    f_q_q(t::Quantity) = sin(0.3∙2π∙t / s) ∙ 9.81N
+    s_1_down  = range(20, 0, length = 100)
+    s_1_up = range(0, 20, length = 100)
+    s_2_down = range(8, 1, length = 100)
+    s_2_up =  range(1, 8, length = 100)
+    plot( xlims =(-5s, 5s),
+        [f_q_q  x-> 0.5f_q_q(1.5x)∙m],
+        ribbon = ([s_1_down∙N  s_2_down∙N∙m], [s_1_up∙N   s_2_up∙N∙m]),
+        yguide = "Load", xguide = "Time")
+end
+
+

test/runtests.jl

@@ -1,6 +1,33 @@
-using MechGluecode
 using Test
+include("plots_1.jl")
+include("diffeq_1.jl")
+include("diffeq_plots_1.jl")
+include("diffeq_plots_2.jl")
 
-@testset "MechGluecode.jl" begin
-    # Write your tests here.
+#=
+using SnoopCompileCore
+list = []
+@testset "Plots" begin
+    push!(list, @snoopr include("plots.jl"))
 end
+@testset "DifferentialEquations" begin
+    push!(list, @snoopr include("differential_equations.jl"))
+end
+using SnoopCompile
+
+invalidations = list[1]
+length(uinvalidated(invalidations))
+trees = invalidation_trees(invalidations)
+length(trees)
+ftrees = filtermod(MechGluecode, trees)
+ftrees = filtermod(MechGlueDiffEq, trees)
+
+
+[(length(trees[i].backedges), i) for i in 1:37] |> sort!
+
+method_invalidations = trees[34]
+length(method_invalidations.backedges)
+root = method_invalidations.backedges[1]
+show(root; minchildren=0)
+
+=#