|
| 1 | +# DefaultInit for all Sundials integrators - handles ModelingToolkit parameter initialization |
| 2 | +function DiffEqBase.initialize_dae!(integrator::AbstractSundialsIntegrator, |
| 3 | + initializealg::DefaultInit) |
| 4 | + prob = integrator.sol.prob |
| 5 | + if prob.f.initialization_data !== nothing |
| 6 | + DiffEqBase.initialize_dae!(integrator, SciMLBase.OverrideInit()) |
| 7 | + else |
| 8 | + DiffEqBase.initialize_dae!(integrator, SciMLBase.CheckInit()) |
| 9 | + end |
| 10 | +end |
| 11 | + |
| 12 | +function DiffEqBase.initialize_dae!(integrator::IDAIntegrator, |
| 13 | + initializealg::BrownFullBasicInit) |
| 14 | + if integrator.u_modified |
| 15 | + IDAReinit!(integrator) |
| 16 | + end |
| 17 | + integrator.f(integrator.tmp, integrator.du, integrator.u, integrator.p, integrator.t) |
| 18 | + tstart, tend = integrator.sol.prob.tspan |
| 19 | + # Use abstol from algorithm parameter (defaults to 1e-10) |
| 20 | + if any(abs.(integrator.tmp) .>= initializealg.abstol) |
| 21 | + if integrator.sol.prob.differential_vars === nothing |
| 22 | + error("Must supply differential_vars argument to DAEProblem constructor to use IDA initial value solver.") |
| 23 | + end |
| 24 | + # BrownFullBasicInit only modifies algebraic variables |
| 25 | + init_type = IDA_YA_YDP_INIT |
| 26 | + # Use preallocated NVector for differential_vars |
| 27 | + if integrator.diff_vars_nvec !== nothing |
| 28 | + integrator.flag = IDASetId(integrator.mem, integrator.diff_vars_nvec) |
| 29 | + else |
| 30 | + error("differential_vars NVector not preallocated but needed for IDASetId") |
| 31 | + end |
| 32 | + dt = integrator.dt == tstart ? tend : integrator.dt |
| 33 | + integrator.flag = IDACalcIC(integrator.mem, init_type, dt) |
| 34 | + |
| 35 | + # Reflect consistent initial conditions back into the integrator's |
| 36 | + # shadow copy. N.B.: ({du, u}_nvec are aliased to {du, u}). |
| 37 | + IDAGetConsistentIC(integrator.mem, integrator.u_nvec, integrator.du_nvec) |
| 38 | + end |
| 39 | + if integrator.t == tstart && integrator.flag < 0 |
| 40 | + integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, |
| 41 | + ReturnCode.InitialFailure) |
| 42 | + end |
| 43 | +end |
| 44 | + |
| 45 | +function DiffEqBase.initialize_dae!(integrator::IDAIntegrator, |
| 46 | + initializealg::ShampineCollocationInit) |
| 47 | + if integrator.u_modified |
| 48 | + IDAReinit!(integrator) |
| 49 | + end |
| 50 | + integrator.f(integrator.tmp, integrator.du, integrator.u, integrator.p, integrator.t) |
| 51 | + tstart, tend = integrator.sol.prob.tspan |
| 52 | + if any(abs.(integrator.tmp) .>= integrator.opts.reltol) |
| 53 | + # ShampineCollocationInit modifies all variables |
| 54 | + init_type = IDA_Y_INIT |
| 55 | + # Use initdt from algorithm if provided, otherwise fall back to integrator.dt |
| 56 | + dt = if initializealg.initdt !== nothing |
| 57 | + initializealg.initdt |
| 58 | + elseif integrator.dt == tstart |
| 59 | + tend |
| 60 | + else |
| 61 | + integrator.dt |
| 62 | + end |
| 63 | + integrator.flag = IDACalcIC(integrator.mem, init_type, dt) |
| 64 | + |
| 65 | + # Reflect consistent initial conditions back into the integrator's |
| 66 | + # shadow copy. N.B.: ({du, u}_nvec are aliased to {du, u}). |
| 67 | + IDAGetConsistentIC(integrator.mem, integrator.u_nvec, integrator.du_nvec) |
| 68 | + end |
| 69 | + if integrator.t == tstart && integrator.flag < 0 |
| 70 | + integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, |
| 71 | + ReturnCode.InitialFailure) |
| 72 | + end |
| 73 | +end |
| 74 | + |
| 75 | +function DiffEqBase.initialize_dae!(integrator::AbstractSundialsIntegrator, |
| 76 | + initializealg::SciMLBase.CheckInit) |
| 77 | + # Not allowed to be a DAE, so no-op |
| 78 | +end |
| 79 | + |
| 80 | +function DiffEqBase.initialize_dae!(integrator::AbstractSundialsIntegrator, initalg::SciMLBase.NoInit) end |
| 81 | + |
| 82 | +function DiffEqBase.initialize_dae!(integrator::AbstractSundialsIntegrator, initalg::SciMLBase.OverrideInit) |
| 83 | + prob = integrator.sol.prob |
| 84 | + nlsolve_alg = KINSOL() |
| 85 | + u0, p, success = SciMLBase.get_initial_values(prob, integrator, prob.f, initalg, Val(isinplace(prob)); nlsolve_alg, abstol = integrator.opts.abstol, reltol = integrator.opts.reltol) |
| 86 | + |
| 87 | + if isinplace(prob) |
| 88 | + integrator.u .= u0 |
| 89 | + if length(integrator.sol.u) == 1 |
| 90 | + integrator.sol.u[1] .= u0 |
| 91 | + end |
| 92 | + else |
| 93 | + integrator.u = u0 |
| 94 | + if length(integrator.sol.u) == 1 |
| 95 | + integrator.sol.u[1] = u0 |
| 96 | + end |
| 97 | + end |
| 98 | + integrator.p = p |
| 99 | + sol = integrator.sol |
| 100 | + @reset sol.prob.p = integrator.p |
| 101 | + integrator.sol = sol |
| 102 | + |
| 103 | + if success |
| 104 | + integrator.u_modified = true |
| 105 | + else |
| 106 | + integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, ReturnCode.InitialFailure) |
| 107 | + end |
| 108 | +end |
| 109 | + |
| 110 | +# Implementation of CheckInit for IDAIntegrator |
| 111 | +function DiffEqBase.initialize_dae!(integrator::IDAIntegrator, |
| 112 | + initializealg::SciMLBase.CheckInit) |
| 113 | + if integrator.u_modified |
| 114 | + IDAReinit!(integrator) |
| 115 | + end |
| 116 | + |
| 117 | + # Evaluate the DAE residual at the initial conditions |
| 118 | + integrator.f(integrator.tmp, integrator.du, integrator.u, integrator.p, integrator.t) |
| 119 | + |
| 120 | + # Check if residuals are within tolerance |
| 121 | + if any(abs.(integrator.tmp) .>= integrator.opts.abstol) |
| 122 | + error(""" |
| 123 | + DAE initialization failed with CheckInit: Initial conditions do not satisfy the DAE constraints. |
| 124 | +
|
| 125 | + The residual norm is $(maximum(abs.(integrator.tmp))), which exceeds the tolerance $(integrator.opts.abstol). |
| 126 | +
|
| 127 | + Note that the initial conditions include both `du0` (derivatives) and `u0` (states), |
| 128 | + and the choice of derivatives must be compatible with the states. |
| 129 | +
|
| 130 | + To resolve this issue, you have several options: |
| 131 | + 1. Fix your initial conditions (both `du0` and `u0`) to satisfy the DAE constraints |
| 132 | + 2. Use Brown's full basic initialization: initializealg = DiffEqBase.BrownFullBasicInit() |
| 133 | + - Optional: specify tolerance with DiffEqBase.BrownFullBasicInit(abstol=1e-8) |
| 134 | + 3. Use Shampine's collocation initialization: initializealg = DiffEqBase.ShampineCollocationInit() |
| 135 | + - Optional: specify initial dt with DiffEqBase.ShampineCollocationInit(0.001) |
| 136 | + 4. If using ModelingToolkit, use: initializealg = SciMLBase.OverrideInit() |
| 137 | +
|
| 138 | + Example for automatic initialization: |
| 139 | + solve(prob, IDA(); initializealg = DiffEqBase.BrownFullBasicInit()) |
| 140 | + """) |
| 141 | + end |
| 142 | +end |
0 commit comments