GitHub

Advanced Examples

Solving Differential Equations

There has been some attempts at passing units through differential equations solvers, and while FlexUnits can do this it comes with certain drawbacks:

  1. Passing units through ODE solvers does take a performance hit (often 2x the time)
  2. Passing units through ODE solvers isn't rigorously tested by their developers and regressions can happen
  3. If done improperly, it's difficult to troubleshoot errors

It also turns out that passing units through an ODE solver is often unnecessary. These solvers already produce correct results as long as the units are coherent (that is, they don't require scaling between dimensions). Unit/dimension errors often show up in setting up the ODE equations/functions, so it makes sense to restrict unit validation activity to this space.

In the section below, we introduce a strategy that can validate units inside differential equations, the benefits of this approach include:

  1. It solves the system at essentially zero added cost
  2. It introduces no added complexity, using the ODE solvers as the developers intended (and test for )
  3. It validates units in a way that is easy to interpret

This strategy involves defining new objects that behave like a user-defined object with units inside the differential equation system, but behave like numerical vectors inside the differential equation solver. This can be done by defining new objects with clear dimensional representations and subtype it to QuantFieldVector abstract class defined in FlexUnits

The problem definition

Let us consider an application where we are using the drag force equation to model a falling object with velocity v and position h. This equation also requires us to consider the following parameters:

  1. Cd the drag force coefficient
  2. A the reference area
  3. ρ the density of the fluid the object is falling through
  4. m the mass of the object
  5. g gravitational acceleration

First, we import necessary packages and define our object as QuantFieldVectors contain state values with static dimensions attached to them. The @D_str macro makes easy work of assigning static dimensions in an intuitive manner.

using FlexUnits, .UnitRegistry
using OrdinaryDiffEqTsit5
using OrdinaryDiffEqRosenbrock
using StaticArrays
using Plots
using BenchmarkTools
using LinearAlgebra

import FlexUnits: QuantFieldVector, DimsMod
import OrdinaryDiffEqTsit5.SciMLBase.FullSpecialize

@kwdef struct FallingObjectState{T} <: QuantFieldVector{2,T}
    v  :: Quantity{T, D"m/s"}
    h  :: Quantity{T, D"m"}
end

@kwdef struct FallingObjectProps{T} <: QuantFieldVector{5,T}
    Cd :: Quantity{T, D""}
    A  :: Quantity{T, D"m^2"}
    ρ  :: Quantity{T, D"kg/m^3"}
    m  :: Quantity{T, D"kg"}
    g  :: Quantity{T, D"m/s^2"}
end

We can write out the equations as we would normally express them; when returning the value, we need to modify the dimensions by 1/s for the derivative. This is done using the DimsMod{D"1/s"} object constructor which modifies a QuantFieldArray object by a static dimension type D"1/s" to produce a time derivative. The constructor takes two arguments:

  1. The original object (FallingObjectState)
  2. A NamedTuple with fields identical to the original object, and values which are quantities having the modified dimensions
function acceleration(u0::AbstractVector, p::FallingObjectProps, t)
    u = FallingObjectState(u0)

    #Drag force
    fd = -sign(u.v)*0.5*p.ρ*u.v^2*p.Cd*p.A
    
    #Drag force effect on state, DimsMod validates the units of the derivative
    dv = (fd/p.m - p.g)
    dh = u.v

    #Validate the units using DimsMod and return the original type (as required by the ODE solver)
    du = DimsMod{D"1/s"}(FallingObjectState, (v=dv, h=dh))
    return convert(typeof(u0), du)
end

The trick with defining the system this way is that calling FallingObjectState() on a unitless vector will apply units to it, and getting indices from it will give you a unitless scalar in SI base units. 

julia> fos = FallingObjectState([1.0,1.0])
2-element FallingObjectState{Float64} with indices SOneTo(2):
 1
 1

julia> fos.v
1.0 m/s

julia> SVector(fos)
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 1.0

Solving the problem

Defining the problem in this manner will give you internal unit verification of your equations (where you need it most) at zero runtime cost regardless as to whether or not you're using an implicit or explicit solver. You can simply solve the system the way you would normally do it

# =============================================================================================================
println("\nExplicit Solution")
# =============================================================================================================
u0 = FallingObjectState{Float64}(v=0.0u"m/s", h=100u"m")
p  = FallingObjectProps{Float64}(Cd=1.0, A=0.1u"m^2", ρ=1.0u"kg/m^3", m=50u"kg", g=9.81u"m/s^2")
abstol = FallingObjectState{Float64}(v=1e-6u"m/s", h=1e-6u"m")
reltol = SA[1e-6, 1e-6]

tspan = dstrip.((0.0u"min", 0.25u"min")) #Time span must be in seconds, dstrip takes care of this
prob = ODEProblem{false, FullSpecialize}(acceleration, u0, tspan, p, abstol=abstol, reltol=reltol)
sol = solve(prob, Tsit5())
plt = plot(sol.t, [dstrip(u.v) for u in sol.u], label="explicit") #Each element in sol.u is a QuantFieldVector

# =============================================================================================================
println("\nImplicit Solution")
# =============================================================================================================
u0 = FallingObjectState{Float64}(v=0.0u"m/s", h=100u"m")
p  = FallingObjectProps{Float64}(Cd=1.0, A=0.1u"m^2", ρ=1.0u"kg/m^3", m=50u"kg", g=9.81u"m/s^2")

tspan = dstrip.((0.0u"min", 0.25u"min"))
prob = ODEProblem{false, FullSpecialize}(acceleration, u0, tspan, p, abstol=abstol, reltol=reltol)
sol = solve(prob, Rodas5P())
plt = plot!(plt, sol.t, [dstrip(u.v) for u in sol.u], label="implicit") #Each element in sol.u is a QuantFieldVector