Merge pull request #886 from Ickaser/dae_init_error

Make CheckInit more informative

Author: ChrisRackauckas

src/SciMLBase.jl

@@ -835,6 +835,8 @@ export ODEFunction, DiscreteFunction, ImplicitDiscreteFunction, SplitFunction, D
 
 export OptimizationFunction, MultiObjectiveOptimizationFunction
 
+export CheckInit
+
 export EnsembleThreads, EnsembleDistributed, EnsembleSplitThreads, EnsembleSerial
 
 export EnsembleAnalysis, EnsembleSummary

src/initialization.jl

@@ -55,11 +55,27 @@ function get_initial_values end
 struct CheckInitFailureError <: Exception
     normresid::Any
     abstol::Any
+    isdae::Bool
 end
 
 function Base.showerror(io::IO, e::CheckInitFailureError)
     print(io,
-        "CheckInit specified but initialization not satisfied. normresid = $(e.normresid) > abstol = $(e.abstol)")
+        "DAE initialization failed: your u0 did not satisfy the initialization requirements, 
+        normresid = $(e.normresid) > abstol = $(e.abstol)."
+    )
+
+    if isdae
+        print(io, " If you wish for the system to 
+            automatically change the algebraic variables to satisfy the algebraic constraints, 
+            please pass `initializealg = BrownBasicInit()` to solve (this option will require 
+            `using OrdinaryDiffEqNonlinearSolve`). If you wish to perform an initialization on the
+            complete u0, please pass initializealg = ShampineCollocationInit() to solve. Note that 
+            initialization can be a very difficult process for DAEs and in many cases can be 
+            numerically intractable without symbolic manipulation of the system. For an automated 
+            system that will generate numerically stable initializations, see ModelingToolkit.jl 
+            structural simplification for more details."
+        )
+    end
 end
 
 struct OverrideInitMissingAlgorithm <: Exception end
@@ -134,7 +150,7 @@ function get_initial_values(
     normresid = isdefined(integrator.opts, :internalnorm) ?
                 integrator.opts.internalnorm(tmp, t) : norm(tmp)
     if normresid > abstol
-        throw(CheckInitFailureError(normresid, abstol))
+        throw(CheckInitFailureError(normresid, abstol, true))
     end
     return u0, p, true
 end
@@ -151,7 +167,7 @@ function get_initial_values(
                 integrator.opts.internalnorm(resid, t) : norm(resid)
 
     if normresid > abstol
-        throw(CheckInitFailureError(normresid, abstol))
+        throw(CheckInitFailureError(normresid, abstol, false))
     end
     return u0, p, true
 end