Add tests

Author: ChrisRackauckas

lib/OrdinaryDiffEqCore/src/initialize_dae.jl

@@ -178,6 +178,14 @@ function _initialize_dae!(integrator, prob::Union{ODEProblem, DAEProblem},
 end
 
 ## CheckInit
+struct CheckInitFailureError <: Exception
+    normresid
+    abstol
+end
+
+function Base.showerror(io::IO, e::CheckInitFailureError)
+    print(io, "CheckInit specified but initialization not satisifed. normresid = $(e.normresid) > abstol = $(e.abstol)")
+end
 
 function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit,
         isinplace::Val{true})
@@ -195,7 +203,7 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit,
 
     normresid = integrator.opts.internalnorm(tmp, t)
     if normresid > integrator.opts.abstol 
-        error("CheckInit specified but initialization not satisifed. normresid = $normresid > abstol = $(integrator.opts.abstol )")
+        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
     end
 end
 
@@ -212,9 +220,9 @@ function _initialize_dae!(integrator, prob::ODEProblem, alg::CheckInit,
     du = f(u0, p, t)
     resid = _vec(du)[algebraic_eqs]
 
-    normresid = integrator.opts.internalnorm(tmp, t)
+    normresid = integrator.opts.internalnorm(resid, t)
     if normresid > integrator.opts.abstol 
-        error("CheckInit specified but initialization not satisifed. normresid = $normresid > abstol = $(integrator.opts.abstol )")
+        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
     end
 end
 
@@ -227,7 +235,7 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     f(resid, integrator.du, u0, p, t)
     normresid = integrator.opts.internalnorm(resid, t)
     if normresid > integrator.opts.abstol 
-        error("CheckInit specified but initialization not satisifed. normresid = $normresid > abstol = $(integrator.opts.abstol )")
+        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
     end
 end
 
@@ -245,6 +253,6 @@ function _initialize_dae!(integrator, prob::DAEProblem,
     resid = f(integrator.du, u0, p, t)
     normresid = integrator.opts.internalnorm(resid, t)
     if normresid > integrator.opts.abstol 
-        error("CheckInit specified but initialization not satisifed. normresid = $normresid > abstol = $(integrator.opts.abstol )")
+        throw(CheckInitFailureError(normresid, integrator.opts.abstol))
     end
 end

test/interface/checkinit_tests.jl

@@ -0,0 +1,50 @@
+using OrdinaryDiffEqBDF, OrdinaryDiffEqRosenbrock, LinearAlgebra, ForwardDiff, Test
+using OrdinaryDiffEqCore
+
+function rober(du, u, p, t)
+    y₁, y₂, y₃ = u
+    k₁, k₂, k₃ = p
+    du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
+    du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
+    du[3] = y₁ + y₂ + y₃ - 1
+    nothing
+end
+function rober(u, p, t)
+    y₁, y₂, y₃ = u
+    k₁, k₂, k₃ = p
+    [-k₁ * y₁ + k₃ * y₂ * y₃,
+        k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2,
+        y₁ + y₂ + y₃ - 1]
+end
+M = [1.0 0 0
+     0 1.0 0
+     0 0 0]
+roberf = ODEFunction(rober, mass_matrix = M)
+roberf_oop = ODEFunction{false}(rober, mass_matrix = M)
+prob_mm = ODEProblem(roberf, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4))
+prob_mm_oop = ODEProblem(roberf_oop, [1.0, 0.0, 0.2], (0.0, 1e5), (0.04, 3e7, 1e4))
+
+@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob_mm, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit())
+@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob_mm_oop, Rodas5P(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit())
+
+f_oop = function (du, u, p, t)
+    out1 = -0.04u[1] + 1e4 * u[2] * u[3] - du[1]
+    out2 = +0.04u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
+    out3 = u[1] + u[2] + u[3] - 1.0
+    [out1, out2, out3]
+end
+
+f = function (resid, du, u, p, t)
+    resid[1] = -0.04u[1] + 1e4 * u[2] * u[3] - du[1]
+    resid[2] = +0.04u[1] - 3e7 * u[2]^2 - 1e4 * u[2] * u[3] - du[2]
+    resid[3] = u[1] + u[2] + u[3] - 1.0
+end
+
+u₀ = [1.0, 0, 0.2]
+du₀ = [0.0, 0.0, 0.0]
+tspan = (0.0, 100000.0)
+differential_vars = [true, true, false]
+prob = DAEProblem(f, du₀, u₀, tspan, differential_vars = differential_vars)
+prob_oop = DAEProblem(f_oop, du₀, u₀, tspan, differential_vars = differential_vars)
+@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit())
+@test_throws OrdinaryDiffEqCore.CheckInitFailureError solve(prob_oop, DFBDF(), reltol = 1e-8, abstol = 1e-8, initializealg = SciMLBase.CheckInit())

test/runtests.jl

@@ -65,7 +65,8 @@ end
         @time @safetestset "Linear Solver Split ODE Tests" include("interface/linear_solver_split_ode_test.jl")
         @time @safetestset "Sparse Diff Tests" include("interface/sparsediff_tests.jl")
         @time @safetestset "Enum Tests" include("interface/enums.jl")
-        @time @safetestset "Enum Tests" include("interface/get_du.jl")
+        @time @safetestset "CheckInit Tests" include("interface/checkinit_tests.jl")
+        @time @safetestset "Get du Tests" include("interface/get_du.jl")
         @time @safetestset "Mass Matrix Tests" include("interface/mass_matrix_tests.jl")
         @time @safetestset "W-Operator prototype tests" include("interface/wprototype_tests.jl")
     end