Skip to content

Commit 245fbbf

Browse files
hustfhustf
committed
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
1 parent 098db0f commit 245fbbf

File tree

9 files changed

+326
-22
lines changed

9 files changed

+326
-22
lines changed

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1 +1,2 @@
11
/Manifest.toml
2+
/.vscode

‎Project.toml‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "MechGluecode"
22
uuid = "3017d99d-ab52-4519-99a2-fa9ddc4637fe"
33
authors = ["hustf <[email protected]> and contributors"]
4-
version = "0.1.4"
4+
version = "0.1.5"
55

66
[deps]
77
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -24,6 +24,8 @@ julia = "1"
2424

2525
[extras]
2626
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
27+
SnoopCompile = "aa65fe97-06da-5843-b5b1-d5d13cad87d2"
28+
SnoopCompileCore = "e2b509da-e806-4183-be48-004708413034"
2729

2830
[targets]
2931
test = ["Test"]

‎README.md‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1 +1,13 @@
11
# MechGluecode
2+
3+
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).
4+
5+
The extension glue code is loaded if you also load:
6+
- DifferentialEquations.jl => load MechGlueDiffEqBase.jl
7+
- Plots.jl => load MechGluePlots.jl
8+
9+
Not completed yet:
10+
- Interpolations.jl => load MechGlueDiffEqBase.jl
11+
- Plots.jl => load MechGlueInterpolations.jl
12+
- ModelingToolkit => MechGlueModelingToolkit
13+
- [RecursiveArrayTools](https://github.com/SciML/RecursiveArrayTools.jl): This is included in MechGlueDiffEqBase, no standalone glue code needed (?)

‎src/MechGluecode.jl‎

Lines changed: 18 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1,37 +1,37 @@
11
module MechGluecode
22
using Requires
3-
export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu
3+
export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu, norm
44

55
function __init__()
66
@require MechanicalUnits = "e6be9192-89dc-11e9-36e6-5dbcb28f419e" begin
77
@require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
88
@info "Plots => using MechGluePlots"
99
@eval using MechGluePlots
1010
end
11-
@require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" begin
12-
import MechanicalUnits: Unitfu, @import_expand
13-
import Unitfu: AbstractQuantity
14-
import DifferentialEquations
15-
import DifferentialEquations: DiffEqBase
16-
import DiffEqBase: value, ODE_DEFAULT_NORM, UNITLESS_ABS2
17-
@info "DifferentialEquations => using MechGlueDiffEqBase"
18-
@eval using MechGlueDiffEqBase
19-
end
20-
2111
@require Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" begin
22-
@info "Interpolations => using MechGlueInterpolations"
23-
@eval using MechGlueInterpolations
12+
#@info "Interpolations => using MechGlueInterpolations"
13+
#@eval using MechGlueInterpolations
2414
end
25-
15+
2616
@require ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" begin
27-
@info "ModelingToolkit => using MechGlueModelingToolkit"
28-
@eval using MechGlueModelingToolkit
17+
#@info "ModelingToolkit => using MechGlueModelingToolkit"
18+
#@eval using MechGlueModelingToolkit
2919
end
3020

3121
@require RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" begin
32-
@eval using MechGlueDiffEqBase
3322
@info "RecursiveArrayTools => using MechGlueRecursiveArrayTools"
34-
@eval using MechGlueRecursiveArrayTools
23+
# @eval using MechGlueDiffEqBase
24+
# @info "RecursiveArrayTools => using MechGlueRecursiveArrayTools"
25+
# @eval using MechGlueRecursiveArrayTools
26+
end
27+
@require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" begin
28+
@info "DifferentialEquations => using MechGlueDiffEqBase"
29+
# @eval import MechanicalUnits: Unitfu, @import_expand
30+
# @eval import Unitfu: AbstractQuantity
31+
# @eval import DifferentialEquations
32+
# @eval import DifferentialEquations: DiffEqBase
33+
# @eval import SciMLBase
34+
@eval using MechGlueDiffEqBase
3535
end
3636
end
3737
@info "MechGluecode init"

‎test/diffeq_1.jl‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
using Test
2+
using MechGluecode
3+
using MechanicalUnits
4+
using MechanicalUnits: Unitfu.promote_to_derived
5+
@import_expand(m, s, N) # Takes long if defined late)
6+
using DifferentialEquations
7+
8+
9+
@testset "Unitful time with unitless state" begin
10+
u0 = 30.0s
11+
tspan = (0.0, 10.0)s
12+
prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
13+
prob2 = ODEProblem((u,t,p) -> (-0.2∙s⁻¹ * u[1]), u0, tspan)
14+
prob3 = ODEProblem((u,t,p) -> [-0.2∙s⁻¹* u[1]], [u0],tspan)
15+
for prob in [prob1, prob2, prob3]
16+
@test solve(prob, Tsit5()).retcode === :Success
17+
end
18+
end
19+
@testset "ODE quantity" begin
20+
f = (y, p, t) -> 0.5y / 3.0s # The derivative, yeah
21+
u0 = 1.5N
22+
tspan = (0.0s, 1.0s)
23+
prob = ODEProblem(f, u0, tspan)
24+
sol = solve(prob)
25+
@test sol.t isa Vector{<:Time}
26+
@test sol.u isa Vector{<:Force}
27+
end
28+

‎test/diffeq_plots_1.jl‎

Lines changed: 120 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,120 @@
1+
using Test
2+
using MechGluecode
3+
using MechanicalUnits
4+
using MechanicalUnits: Unitfu.promote_to_derived
5+
@import_expand(m, s, N) # Takes long if defined late)
6+
using DifferentialEquations
7+
8+
9+
@testset "Unitful time with unitful function" begin
10+
u0 = 30.0N
11+
tspan = (0.0, 10.0)s
12+
prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
13+
prob2 = ODEProblem((u,t,p) -> (-0.2∙s⁻¹ * u[1]), u0, tspan)
14+
prob3 = ODEProblem((u,t,p) -> [-0.2∙s⁻¹* u[1]], [u0],tspan)
15+
so1 = solve(prob1, Euler(), dt = 1s)
16+
so2 = solve(prob2, Midpoint())
17+
so3 = solve(prob3, Tsit5())
18+
plot(so1, denseplot = false, marker = true)
19+
plot(so2, denseplot = false, marker = true)
20+
plot(so3, denseplot = false, marker = true)
21+
sf1(t) = so1(t)[1]
22+
sf2(t) = so2(t)[1]
23+
sf3(t) = so3(t)[1]
24+
plot([sf1, sf2, sf3], xlims =(0,10)s, ylims =(0,40)N)
25+
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))
26+
27+
end
28+
@testset "Unitful time with unitless function" begin
29+
u0 = 30.0
30+
tspan = (0.0, 10.0)s
31+
prob1 = ODEProblem((du,u,t,p) -> (du[1] = -0.2∙s⁻¹ * u[1]), [u0], tspan)
32+
prob2 = ODEProblem((u,t,p) -> (-0.2∙s⁻¹ * u[1]), u0, tspan)
33+
prob3 = ODEProblem((u,t,p) -> [-0.2∙s⁻¹* u[1]], [u0],tspan)
34+
so1 = solve(prob1, Euler(), dt = 1s)
35+
so2 = solve(prob2, Midpoint())
36+
so3 = solve(prob3, Tsit5())
37+
# don't work plot(so1, denseplot = false, marker = true)
38+
# don't work plot(so2, denseplot = false, marker = true)
39+
# don't work plot(so3, denseplot = false, marker = true)
40+
sf1(t) = so1(t)[1]
41+
sf2(t) = so2(t)[1]
42+
sf3(t) = so3(t)[1]
43+
plot([sf1, sf2, sf3], xlims =(0,10)s, ylims =(0,40))
44+
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))
45+
end
46+
47+
@testset "Plot ODE quantity" begin
48+
# Don't give us kN from from solution in this case, just defaults.
49+
promote_to_derived()
50+
f1 = (y, p, t) -> -1.5y / 0.3s
51+
u0 = 1.5N
52+
tspan = (0.0s, 1.0s)
53+
prob = ODEProblem(f1, u0, tspan)
54+
sol = solve(prob)
55+
plot(sol)
56+
end
57+
58+
@testset "Plot ODE system same quantities" begin
59+
# This is not a good example, as the units are incorrect and have no
60+
# known physical interpretation
61+
σ = 10/s
62+
ρ = 28m
63+
β = (8 / 3)m
64+
function lorenz!(du, u, p, t)
65+
du[1] = (u[2] - u[1])∙σ
66+
du[2] = (u[1] * (ρ - u[3])) / (m∙s) - u[2]/s
67+
du[3] = u[1] * u[2] / (m∙s) - β ∙ u[3]/(m∙s)
68+
end
69+
u0 = [1.0, 0.0, 0.0]m
70+
tspan = (0.0,100.0)s
71+
prob = ODEProblem(lorenz!,u0,tspan)
72+
sol = solve(prob)
73+
# Plot all variables vs time
74+
plot(sol)
75+
# Plot in phase space
76+
plot(sol, vars=(1,2,3))
77+
end
78+
#=
79+
@testset "Plot ODE system single pendulum" begin
80+
# Second order non-linear, undamped pendulum:
81+
# mLθ'' + mg∙sin(θ)
82+
# where m is mass, unit kg
83+
# L is length, unit m
84+
# θ is angle to vertical, strictly no unit
85+
# g is gravity's acceleration
86+
# Acceleration
87+
# θ'' + k * sin(θ) = 0
88+
# where k = g / L , unit(k) = s⁻²
89+
# Reduced to first order non-linear system by defining:
90+
# θ = y[1]
91+
# θ' = y[2]
92+
# It follows that
93+
# θ'' = y'[2]
94+
# and
95+
# y'[1] = y[2]
96+
# y'[2] = -k * sin(y[1])
97+
# To solve numerically, we define the derivative.
98+
# For efficiency, it's evaluated in place, i.e. dy changes at every call.
99+
function f!(dy, y, k, t)
100+
dy[1] = y[2]
101+
dy[2] = -k ∙ sin(y[1])
102+
end
103+
# Fixed parameter
104+
k = upreferred(g / 1.0m)
105+
# Initial conditions
106+
θ0 = 0.1
107+
θ0_vel = 0.0s⁻¹
108+
θ0_acc = -k *sin(θ0)
109+
110+
y0 = [θ0, θ0_vel]
111+
y0 = ArrayPartition([θ0_vel, θ0_acc])
112+
tspan = (0.0, 100.0)s
113+
prob = ODEProblem(f!, y0, tspan, k)
114+
sol = solve(prob)
115+
# Plot all variables vs time
116+
plot(sol)
117+
# Plot in phase space
118+
plot(solve(prob), vars=(1,2))
119+
end
120+
=#

‎test/diffeq_plots_2.jl‎

Lines changed: 84 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,84 @@
1+
using Test
2+
using MechGluecode
3+
using MechanicalUnits
4+
using MechanicalUnits: Unitfu.promote_to_derived
5+
@import_expand(~m, s, N)
6+
using DifferentialEquations
7+
import MechGluecode: norm # should work without
8+
9+
@testset "Initial checks ArrayPartition" begin
10+
r0 = [1131.340, -2282.343, 6672.423]∙km
11+
v0 = [-5.64305, 4.30333, 2.42879]∙km/s
12+
Δt = 86400.0*365∙s
13+
μ = 398600.4418∙km³/s²
14+
rv0 = ArrayPartition(r0,v0)
15+
16+
function fo(dy, y, μ, t)
17+
r = norm(y.x[1])
18+
dy.x[1] .= y.x[2]
19+
dy.x[2] .= -μ .* y.x[1] / r^3
20+
end
21+
prob = ODEProblem(fo, rv0, (0.0s, 0.035Δt), μ)
22+
@time sol = solve(prob, Tsit5(), dt = 0.1s);
23+
# Plot all variables vs time
24+
plot(sol)
25+
plot(sol, vars = (1))
26+
plot(sol, vars = (2))
27+
plot(sol, vars = (3))
28+
plot(sol, vars = (4))
29+
plot(sol, vars = (5))
30+
plot(sol, vars = (6))
31+
32+
# Plot in phase space
33+
plot(sol, vars=(1,4))
34+
plot(sol, vars=(2,5))
35+
plot(sol, vars=(3,6))
36+
plot(sol, vars=(2,5))
37+
# Three distances vs time (0)
38+
plot(sol, vars=[(0,1), (0,2),(0,3)])
39+
# Three velocities vs time (0)
40+
plot(sol, vars=[(0,4), (0,5),(0,6)])
41+
end
42+
#=
43+
@testset "Plot ODE system single pendulum" begin
44+
# Second order non-linear, undamped pendulum:
45+
# mLθ'' + mg∙sin(θ)
46+
# where m is mass, unit kg
47+
# L is length, unit m
48+
# θ is angle to vertical, strictly no unit
49+
# g is gravity's acceleration
50+
# Acceleration
51+
# θ'' + k * sin(θ) = 0
52+
# where k = g / L , unit(k) = s⁻²
53+
# Reduced to first order non-linear system by defining:
54+
# θ = y[1]
55+
# θ' = y[2]
56+
# It follows that
57+
# θ'' = y'[2]
58+
# and
59+
# y'[1] = y[2]
60+
# y'[2] = -k * sin(y[1])
61+
# To solve numerically, we define the derivative.
62+
# For efficiency, it's evaluated in place, i.e. dy changes at every call.
63+
function f!(dy, y, k, t)
64+
dy[1] = y[2]
65+
dy[2] = -k ∙ sin(y[1])
66+
end
67+
# Fixed parameter
68+
k = upreferred(g / 1.0m)
69+
# Initial conditions
70+
θ0 = 0.1
71+
θ0_vel = 0.0s⁻¹
72+
θ0_acc = -k *sin(θ0)
73+
74+
y0 = [θ0, θ0_vel]
75+
y0 = ArrayPartition([θ0_vel, θ0_acc])
76+
tspan = (0.0, 100.0)s
77+
prob = ODEProblem(f!, y0, tspan, k)
78+
sol = solve(prob)
79+
# Plot all variables vs time
80+
plot(sol)
81+
# Plot in phase space
82+
plot(solve(prob), vars=(1,2))
83+
end
84+
=#

‎test/plots_1.jl‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
using Test
2+
import MechanicalUnits: @import_expand
3+
@import_expand(m, s, N) # Takes long if defined late)
4+
using MechGluecode
5+
import Plots
6+
import MechGluecode
7+
8+
@testset "Plots mixed quantity" begin
9+
s1x = [1.0,2,3]
10+
s2x = [1.0,2,3.5]
11+
s1y = [1.0,2,4]
12+
s2y = [1.0,2,2]
13+
mxqm = hcat(s1x∙m, s2x∙m²)
14+
myqm = hcat(s1y∙s, s2y∙s²)
15+
plot(mxqm, myqm, seriestype = [:path :scatter], label = ["path" "scatter"], xguide = "x", yguide = "y")
16+
end
17+
18+
@testset "Plots functions with mixed units" begin
19+
f_q_q(t::Quantity) = sin(0.3∙2π∙t / s) ∙ 9.81N
20+
s_1_down = range(20, 0, length = 100)
21+
s_1_up = range(0, 20, length = 100)
22+
s_2_down = range(8, 1, length = 100)
23+
s_2_up = range(1, 8, length = 100)
24+
plot( xlims =(-5s, 5s),
25+
[f_q_q x-> 0.5f_q_q(1.5x)∙m],
26+
ribbon = ([s_1_down∙N s_2_down∙N∙m], [s_1_up∙N s_2_up∙N∙m]),
27+
yguide = "Load", xguide = "Time")
28+
end
29+
30+

‎test/runtests.jl‎

Lines changed: 30 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,33 @@
1-
using MechGluecode
21
using Test
2+
include("plots_1.jl")
3+
include("diffeq_1.jl")
4+
include("diffeq_plots_1.jl")
5+
include("diffeq_plots_2.jl")
36

4-
@testset "MechGluecode.jl" begin
5-
# Write your tests here.
7+
#=
8+
using SnoopCompileCore
9+
list = []
10+
@testset "Plots" begin
11+
push!(list, @snoopr include("plots.jl"))
612
end
13+
@testset "DifferentialEquations" begin
14+
push!(list, @snoopr include("differential_equations.jl"))
15+
end
16+
using SnoopCompile
17+
18+
invalidations = list[1]
19+
length(uinvalidated(invalidations))
20+
trees = invalidation_trees(invalidations)
21+
length(trees)
22+
ftrees = filtermod(MechGluecode, trees)
23+
ftrees = filtermod(MechGlueDiffEq, trees)
24+
25+
26+
[(length(trees[i].backedges), i) for i in 1:37] |> sort!
27+
28+
method_invalidations = trees[34]
29+
length(method_invalidations.backedges)
30+
root = method_invalidations.backedges[1]
31+
show(root; minchildren=0)
32+
33+
=#

0 commit comments

Comments
 (0)