Don't check units during initialization

...because initialization problems don't correctly handle units and units have already been checked during system initialization.

Author: ctessum

src/systems/problem_utils.jl

@@ -396,7 +396,7 @@ Keyword arguments:
 - `eval_module`: If `eval_expression == true`, the module to `eval` into. Otherwise, the module
   in which to generate the `RuntimeGeneratedFunction`.
 - `fully_determined`: Override whether the initialization system is fully determined.
-- `check_units`: Enable or disable unit checks.
+- `check_units`: Enable or disable unit checks when constructing the initialization problem.
 - `tofloat`, `use_union`: Passed to [`better_varmap_to_vars`](@ref) for building `u0` (and
   possibly `p`).
 - `u0_constructor`: A function to apply to the `u0` value returned from `better_varmap_to_vars`
@@ -414,7 +414,7 @@ function process_SciMLProblem(
         implicit_dae = false, t = nothing, guesses = AnyDict(),
         warn_initialize_determined = true, initialization_eqs = [],
         eval_expression = false, eval_module = @__MODULE__, fully_determined = false,
-        check_units = true, tofloat = true, use_union = false,
+        check_units = false, tofloat = true, use_union = false,
         u0_constructor = identity, du0map = nothing, check_length = true, symbolic_u0 = false, kwargs...)
     dvs = unknowns(sys)
     ps = parameters(sys)

test/initializationsystem.jl

@@ -3,6 +3,7 @@ using ForwardDiff
 using SymbolicIndexingInterface, SciMLStructures
 using SciMLStructures: Tunable
 using ModelingToolkit: t_nounits as t, D_nounits as D
+using DynamicQuantities
 
 @parameters g
 @variables x(t) y(t) [state_priority = 10] λ(t)
@@ -860,3 +861,19 @@ end
     integ = init(prob)
     @test integ[x] ≈ [1.0, 3.0]
 end
+
+@testset "units" begin
+    t = ModelingToolkit.t
+    D = ModelingToolkit.D
+    @parameters g [unit = u"m/s^2"] L=1 [unit = u"m^2"]
+    @variables x(t) [unit = u"m"] y(t) [unit = u"m" state_priority = 10] λ(t) [unit = u"s^-2"]
+    eqs = [D(D(x)) ~ λ * x
+           D(D(y)) ~ λ * y - g
+           x^2 + y^2 ~ L]
+    @mtkbuild pend = ODESystem(eqs, t)
+
+    prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1],
+        guesses = ModelingToolkit.missing_variable_defaults(pend))
+    sol = solve(prob, Rodas5P())
+    @test SciMLBase.successful_retcode(sol)
+end