Skip to content

Commit

Permalink
modified: .gitignore
Browse files Browse the repository at this point in the history 
modified:   Project.toml DiffEqBase, RecursiveArrayTools, Unitfu only
modified:   README.md    Short explanation
modified:   src/MechGlueDiffEqBase.jl Drop experiments, norm only.
deleted:    test/Manifest.toml
deleted:    test/Project.toml
modified:   test/runtests.jl  Include 1..3
new file:   test/test_1.jl    First version
new file:   test/test_2.jl    First version
new file:   test/test_3.jl    First version
  • Loading branch information
@hustf
hustf committed May 1, 2021
1 parent cade4f2 commit 772a0fc
Show file tree
Hide file tree
Showing 10 changed files with 224 additions and 1,535 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
/Manifest.toml
/.vscode
8 changes: 3 additions & 5 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,21 +1,19 @@
name = "MechGlueDiffEqBase"
uuid = "2532746b-52b5-4539-9431-8bb183ab067f"
authors = ["hustf <[email protected]> and contributors"]
version = "0.1.31"
version = "0.1.32"

[deps]
DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
MechGluePlots = "03a24088-0a82-4850-846f-2cad5d056ae1"
MechanicalUnits = "e6be9192-89dc-11e9-36e6-5dbcb28f419e"
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
Unitfu = "5ee08b94-2369-4f4a-b8c7-99333ba35fb0"

[compat]
julia = "1"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
MechanicalUnits = "e6be9192-89dc-11e9-36e6-5dbcb28f419e"

[targets]
test = ["Test"]
14 changes: 14 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1 +1,15 @@
# MechGlueDiffEqBase
Glue code for making [DiffEqBase](https://github.com/SciML/DiffEqBase.jl) work with units.

It also includes glue code for [RecursiveArrayTools](https://github.com/SciML/RecursiveArrayTools.jl), which enables type-stable solution of equations with mixed units.

This defines how to calculate the vector norm when the vector is given in units compatible with [Unitfu.jl](https://github.com/hustf/Unitfu.jl), from registy [M8](https://github.com/hustf/M8). The differential equation algorithms expects the norm to be unitless, as can be seen in e.g. step size estimators:

err_scaled = **error** / (**abstol** + norm(u) * **reltol**)

where **bold** indicates unitful objects.

The functions are adaptions of corresponding code from [DiffEqBase](https://github.com/SciML/DiffEqBase.jl/blob/6bb8830711e729ef513f2b1beb95853e4a691375/src/init.jl).



144 changes: 21 additions & 123 deletions src/MechGlueDiffEqBase.jl
View file
Original file line number Diff line number Diff line change
@@ -1,139 +1,37 @@
module MechGlueDiffEqBase
import Unitfu: AbstractQuantity, Quantity, ustrip
import DiffEqBase: value, ODE_DEFAULT_NORM, UNITLESS_ABS2
import Unitfu: AbstractQuantity, Quantity, ustrip, norm
import DiffEqBase: value, ODE_DEFAULT_NORM, UNITLESS_ABS2, zero
import DiffEqBase: calculate_residuals, @muladd
using RecursiveArrayTools
export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu, AbstractQuantity, Quantity
export norm, ArrayPartition


import RecipesBase
import RecipesBase.@recipe
import SciMLBase
import SciMLBase: AbstractTimeseriesSolution, RecipesBase
import SciMLBase: AbstractDiscreteProblem, AbstractRODESolution, SensitivityInterpolation
import SciMLBase: getsyms, interpret_vars, cleansyms
import SciMLBase: diffeq_to_arrays, issymbollike, getindepsym_defaultt
import SciMLBase: DEFAULT_PLOT_FUNC


# This is identical to what DiffEqBase defines for Unitful
function value(x::Type{AbstractQuantity{T,D,U}}) where {T,D,U}
T
end
function value(x::AbstractQuantity)
x.val
end
# This is different from what DiffEqBase defines for Unitful
value(::Type{<:AbstractQuantity{T,D,U}}) where {T,D,U<:Core.TypeofBottom} = Base.undef_ref_str
value(x::Q) where {Q<:AbstractQuantity} = ustrip(x)


# This is identical to what DiffEqBase defines for Unitful
@inline function ODE_DEFAULT_NORM(u::AbstractArray{<:AbstractQuantity,N},t) where {N}
# Support adaptive errors should be errorless for exponentiation
sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
# This is identical to what DiffEqBase defines for Unitful
@inline function ODE_DEFAULT_NORM(u::Array{<:AbstractQuantity,N},t) where {N}
sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::AbstractQuantity,t)
abs(value(u))

# This is identical to what DiffEqBase defines for Unitful

@inline function ODE_DEFAULT_NORM(u::AbstractQuantity, t)
abs(ustrip(u))
end
# This is slightly different from what DiffEqBase defines for Unitful
@inline function UNITLESS_ABS2(x::AbstractQuantity)
real(abs2(x)/oneunit(x)*oneunit(x))
real(abs2(x) / (oneunit(x)^2))
end

# This recipe is copied from SciMLBase, and only the signature is adopted
# for use with unitful solutions. It can perhaps be deleted now, and the link to RecipesBase be deleted.
@recipe function f(sol::AbstractTimeseriesSolution{T, N, A};
plot_analytic=false,
denseplot = (sol.dense ||
typeof(sol.prob) <: AbstractDiscreteProblem) &&
!(typeof(sol) <: AbstractRODESolution) &&
!(hasfield(typeof(sol),:interp) &&
typeof(sol.interp) <: SensitivityInterpolation),
plotdensity = min(Int(1e5),sol.tslocation==0 ?
(typeof(sol.prob) <: AbstractDiscreteProblem ?
max(1000,100*length(sol)) :
max(1000,10*length(sol))) :
1000*sol.tslocation),
tspan = nothing, axis_safety = 0.1,
vars=nothing) where {T<:Quantity, N, A<:Array{<:Quantity}}
syms = getsyms(sol)
int_vars = interpret_vars(vars,sol,syms)
strs = cleansyms(syms)

tscale = get(plotattributes, :xscale, :identity)
plot_vecs,labels = diffeq_to_arrays(sol,plot_analytic,denseplot,
plotdensity,tspan,axis_safety,
vars,int_vars,tscale,strs)

tdir = sign(sol.t[end]-sol.t[1])
xflip --> tdir < 0
seriestype --> :path

# Special case labels when vars = (:x,:y,:z) or (:x) or [:x,:y] ...
if typeof(vars) <: Tuple && (issymbollike(vars[1]) && issymbollike(vars[2]))
xguide --> issymbollike(int_vars[1][2]) ? Symbol(int_vars[1][2]) : strs[int_vars[1][2]]
yguide --> issymbollike(int_vars[1][3]) ? Symbol(int_vars[1][3]) : strs[int_vars[1][3]]
if length(vars) > 2
zguide --> issymbollike(int_vars[1][4]) ? Symbol(int_vars[1][4]) : strs[int_vars[1][4]]
end
end

if (!any(issymbollike,getindex.(int_vars,1)) && getindex.(int_vars,1) == zeros(length(int_vars))) ||
(!any(issymbollike,getindex.(int_vars,2)) && getindex.(int_vars,2) == zeros(length(int_vars))) ||
all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,1)) || all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,2))
xguide --> "$(getindepsym_defaultt(sol))"
end
if length(int_vars[1]) >= 3 && ((!any(issymbollike,getindex.(int_vars,3)) && getindex.(int_vars,3) == zeros(length(int_vars))) ||
all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,3)))
yguide --> "$(getindepsym_defaultt(sol))"
end
if length(int_vars[1]) >= 4 && ((!any(issymbollike,getindex.(int_vars,4)) && getindex.(int_vars,4) == zeros(length(int_vars))) ||
all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,4)))
zguide --> "$(getindepsym_defaultt(sol))"
end

if (!any(issymbollike,getindex.(int_vars,2)) && getindex.(int_vars,2) == zeros(length(int_vars))) ||
all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,2))
if tspan === nothing
if tdir > 0
xlims --> (sol.t[1],sol.t[end])
else
xlims --> (sol.t[end],sol.t[1])
end
else
xlims --> (tspan[1],tspan[end])
end
else
mins = minimum(sol[int_vars[1][2],:])
maxs = maximum(sol[int_vars[1][2],:])
for iv in int_vars
mins = min(mins,minimum(sol[iv[2],:]))
maxs = max(maxs,maximum(sol[iv[2],:]))
end
xlims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)
end

# Analytical solutions do not save enough information to have a good idea
# of the axis ahead of time
# Only set axis for animations
if sol.tslocation != 0 && !(typeof(sol) <: AbstractAnalyticalSolution)
if all(getindex.(int_vars,1) .== DEFAULT_PLOT_FUNC)
mins = minimum(sol[int_vars[1][3],:])
maxs = maximum(sol[int_vars[1][3],:])
for iv in int_vars
mins = min(mins,minimum(sol[iv[3],:]))
maxs = max(maxs,maximum(sol[iv[3],:]))
end
ylims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)

if length(int_vars[1]) >= 4
mins = minimum(sol[int_vars[1][4],:])
maxs = maximum(sol[int_vars[1][4],:])
for iv in int_vars
mins = min(mins,minimum(sol[iv[4],:]))
maxs = max(mins,maximum(sol[iv[4],:]))
end
zlims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)
end
end
end

label --> reshape(labels,1,length(labels))
println("----------------------------------MechGlueDiffEqBase")
(plot_vecs...,)
end

end
Loading

0 comments on commit 772a0fc

Please sign in to comment.