modified: README.md Point out test_5

modified:   src/MechGlueDiffEqBase.jl  Two-param 'similar'
modified:   test/runtests.jl       test_5.jl
modified:   test/test_2.jl         spacing
modified:   test/test_3.jl         Comment inferrability
modified:   test/test_4.jl         Test two-param 'similar'
new file:   test/test_5.jl         Inferrable mixed-units calc

Author: hustf

README.md

@@ -11,6 +11,7 @@ These functions preserve types for mixed-unit vectors. E.g. [1.0kg, 2.0N, 3.0m/s
 
 This package also uses and reexports 'ArrayPartition' from [RecursiveArrayTools](https://github.com/SciML/RecursiveArrayTools.jl), which enables type-stable solution of equations with mixed units. It pre-compiles it with use of mixed unit vectors.
 
+An example of inferrable, mixed units calculations with debugging is included in `test/test_5.jl`.
 
 Note: The way error tolerances are defined is initially confusing, but good to know:
 

src/MechGlueDiffEqBase.jl

@@ -61,16 +61,23 @@ end
 
 # Vectors with compatible units, treat as normal
 similar(x::Vector{Quantity{T, D, U}}) where {T,D,U} = Vector{Quantity{T, D, U}}(undef, size(x,1))
-#similar(a::Array{T,1}) where {T}                    = Vector{T}(undef, size(a,1))
 # Vectors with incompatible units, special inferreable treatment
 # VERY similar is still similar and (very slightly) different
 similar(x::Vector{Q}) where {Q<:AbstractQuantity{T, D, U} where {D,U}} where T = copy(x)
 
 
+# Vectors with compatible units, treat as normal
+similar(x::Vector{Quantity{T, D, U}}, S::Type) where {T,D,U} = Vector{S}(undef, size(x,1))
+# Vectors with incompatible units, special inferreable treatment
+function similar(x::Vector{Q}, S::Type) where {Q<:AbstractQuantity{T, D, U} where {D, U}} where T
+    x0 = Vector{S}(undef, size(x, 1))
+end
+
 
 
 # KISS pre-compillation to reduce loading times
 # This is simply a boiled-down obfuscated test_4.jl
+
 import Unitfu: m, s, kg, N, ∙
 let
     r0ul = [1131.340, -2282.343, 6672.423]

test/runtests.jl

@@ -2,4 +2,5 @@ using Test
 include("test_1.jl")
 include("test_2.jl")
 include("test_3.jl")
-include("test_4.jl")
+include("test_4.jl")
+include("test_5.jl")

test/test_2.jl

@@ -46,12 +46,6 @@ end
     sol = solve(prob,ExplicitRK())
 end
 
-
-
-
-
-
-
 @testset "Without ArrayPartition - broken" begin
     # coordinate: u = [position, momentum]
     # parameters: p = [mass, force constanst]

test/test_3.jl

@@ -70,7 +70,7 @@ end
 
 
 
-@testset "ArrayPartition with mixed units" begin
+@testset "ArrayPartition with mixed units, not inferrable" begin
     α₀() = 30°
     x₀() = 0.0m
     y₀() = 0.0m
@@ -88,6 +88,8 @@ end
     Rx(vx, vy) = R(vx, vy) * cos(α(vx, vy))
     Ry(vx, vy) = R(vx, vy) * sin(α(vx, vy))
 
+    # An inferrable version of this is much faster
+    # but this should work anyway.
     function f(du,u,p,t)
         x, y, vx,vy = u
         du[1] = dx = vx

test/test_4.jl

@@ -3,6 +3,7 @@ using Test
 using MechGlueDiffEqBase
 using MechanicalUnits: @import_expand, dimension, NoDims, ∙
 import MechanicalUnits: g, g⁻¹
+import MechanicalUnits.Unitfu.numtype
 @import_expand(km, N, s, m, km, kg, °, inch)
 using DiffEqBase, OrdinaryDiffEq
 
@@ -63,6 +64,9 @@ end
     sima = @inferred similar(rv0)
     @test all(typeof.(rv0) == typeof.(sima))
     @test sima !== rv0
+    # The expected behaviour (in some algorithms) is to drop the units here:
+    simb = @inferred(similar(r0, Int64))
+    @test all(typeof.(simb) .== Int64)
 end
 
 
@@ -75,6 +79,8 @@ end
     sima = @inferred similar(rv0)
     @test all(typeof.(rv0) == typeof.(sima))
     @test sima !== rv0
+    simb = @inferred(similar(r0, Int64))
+    @test all(typeof.(simb) .== Int64)
 end
 
 @testset "Inferrable UNITLESS_ABS2 Unitless ArrayPartition" begin

test/test_5.jl

@@ -0,0 +1,80 @@
+# Differential equations with mixed units using ArrayPartition and no specified algorithm
+using MechGlueDiffEqBase
+using MechanicalUnits: @import_expand, ∙, g
+@import_expand(km, N, s, m, km, kg, °, inch)
+using OrdinaryDiffEq: ODEProblem
+using DifferentialEquations: solve, Tsit5
+using Logging: disable_logging, LogLevel, Debug, with_logger
+import Logging
+
+using Test: @inferred
+
+# Constants we don't think we'll change ever
+const x₀ = 0.0km
+const y₀ = 0.0km
+const ø = 15inch
+const ρ = 1.225kg/m³
+
+# Calculated constants
+const Aₚᵣ = π/4 * ø^2
+# Constants that we define as functions, because we may
+# want to modify them later in the same scripting session.
+α₀()  = 30°
+v₀()  = 1050m/s
+
+mₚ()   = 495kg
+C_s() = 0.4
+x´₀() = v₀() * cos(α₀())
+y´₀() = v₀() * sin(α₀())
+
+# Local tuple initial condition
+u₀ = ArrayPartition([x₀, y₀], [x´₀(), y´₀()])
+
+# Functions
+v(vx, vy) = √(vx^2 + vy^2)
+R(vx, vy) = 0.5∙ρ∙C_s()∙Aₚᵣ∙v(vx, vy)^2
+α(vx, vy) = atan(vy, vx)
+Rx(vx, vy) = R(vx, vy) * cos(α(vx, vy))
+Ry(vx, vy) = R(vx, vy) * sin(α(vx, vy))
+
+
+# Local tuple, i.e. the interesting degrees of freedom
+# and their derivatives
+function  Γ!(u´, u, p, t)
+    (x, y), (x´, y´) = packout(u)
+    @debug "Γ"    x    y    x´   y´   maxlog = 3
+    # Calculate the acceleration for this step
+    x´´ =     -Rx(x´, y´) / mₚ()
+    y´´ = -1g -Ry(x´, y´) / mₚ()
+    @debug "Γ Acceleration"    x´´|> g    y´´|>g    maxlog = 3
+    packin!(u´, x´, y´, x´´, y´´)
+end
+
+packout(u) = u.x[1], u.x[2]
+function packin!(u´, x´, y´, x´´, y´´)
+    # We access the values through its field name, x (unrelated to the specific problem)
+    u´.x[1][1] = x´;  u´.x[1][2] = y´
+    u´.x[2][1] = x´´; u´.x[2][2] = y´´
+    return u´
+end
+
+function solve_guarded(u₀; alg = Tsit5(), debug=false)
+    # Test the functions
+    @inferred packout(u₀)
+    @inferred packout(0.01u₀/s)
+    @inferred Γ!(u₀/s,u₀, nothing, nothing)
+    !debug  && disable_logging(Debug)
+    prob = ODEProblem( Γ!,u₀,(0.0,60)s)
+    sol = if debug
+        with_logger(Logging.ConsoleLogger(stderr, Logging.Debug)) do
+            solve(prob, Tsit5())
+        end
+    else
+        solve(prob, Tsit5())
+    end
+    # Re-enable
+    disable_logging(LogLevel(Debug-1))
+    sol
+end
+@time sol = solve_guarded(u₀, debug=true);
+nothing