Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add new option to DAE solvers to never initialize DAE #549

Closed
MartinOtter opened this issue Jan 16, 2020 · 9 comments
Closed

Add new option to DAE solvers to never initialize DAE #549

MartinOtter opened this issue Jan 16, 2020 · 9 comments

Comments

@MartinOtter
Copy link

Please, add a new option to the DAE solvers, so that initialization is never performed by the DAE solver.

Note, the IDA DAE solver uses the following code in Sundials/src/common_interface/solve.jl:

    f!(rtest, du0, u0, prob.p, t0)
    if any(abs.(rtest) .>= reltol)
        ...
        flag = IDACalcIC(mem, init_type, _t)
    end 

This means that solve first calls the user-suppled function f! = prob.f and then checks whether the absolute values of the returned residuals rtest are small enough. If this is not the case, the IDA initialization IDACalcIC is called. The problem is that this function may easily fail, even if the IDA integrator will integrate the model with the provided u0, du0. Explanation:

Initialization of DAEs is difficult. Good initializers use model information that is not available to a general DAE solver. For example assume that the DAE is initialized in steady state (so du=0) with a homotopy method (because good start values are not available for u0; this is, for example, standard for electric circuit simulators). Good non-linear solvers will not use the absolute values of the residuals as termination conditions, because this is not scale invariant (multiplying a residual equation with any number, will not change the exact mathematical solution). Therefore, it is easily possible that the solver test any(abs.(rtest) .>= reltol) fails, even if initialization is fine and the DAE integrator will integrate the model with the u0, du0 computed by an external initializer.

@ChrisRackauckas
Copy link
Member

Indeed I agree this option is necessary.

@ChrisRackauckas
Copy link
Member

We will also want it accessible in callbacks so after modifying values you can tell it whether or not to re-run initialization.

@MartinOtter
Copy link
Author

I could not figure out what to do (with DifferentialEquations v6.15.0). I tried to call

DifferentialEquations.u_modified!(integrator,false) 

before running a simulation and then get the message

ERROR: LoadError: MethodError: no method matching u_modified!(::Sundials.IDA{:Dense,Nothing,Nothing}, ::Bool)
Closest candidates are:
  u_modified!(::OrdinaryDiffEq.ODEIntegrator, ::Bool) at D:\otter\home\.julia\packages\OrdinaryDiffEq\NsugH\src\integrators\integrator_interface.jl:48
  u_modified!(::StochasticDiffEq.SDEIntegrator, ::Bool) at D:\otter\home\.julia\packages\StochasticDiffEq\BhcRo\src\integrators\integrator_interface.jl:32
  u_modified!(::DelayDiffEq.DDEIntegrator, ::Bool) at D:\otter\home\.julia\packages\DelayDiffEq\BCL4Q\src\integrators\interface.jl:157
  ...

It seems that u_modified! is not supported for a DAE integrator?

Since you closed this issue, it seems to be solved. Can you give more explanations what to do to switch off automatic initialization for DAE problems

@ChrisRackauckas
Copy link
Member

Oh, forgot to add it to the docs: https://github.com/SciML/DiffEqDocs.jl/blob/master/docs/src/solvers/dae_solve.md#full-list-of-methodsid-dae_solve_full .

It seems that u_modified! is not supported for a DAE integrator?

It is. https://github.com/SciML/Sundials.jl/blob/master/src/common_interface/integrator_utils.jl#L113-L115 What version of Sundials are you using?

@ChrisRackauckas
Copy link
Member

Now that NoInit() is in the docs I'll close this.

(We still need to unify Sundials with the new initialization setup, but for now it's 🤷 )

@MartinOtter
Copy link
Author

Still do not manage. I tried:

import DifferentialEquations
DifferentialEquations.solve(<..>, initializealg = DifferentialEquations.NoInit())   

and get then the error

ERROR: LoadError: UndefVarError: NoInit not defined

I also tried to call u_modified! directly before calling solve

import DifferentialEquations
integrator = DifferentialEquations.IDA()
DifferentialEquations.u_modified!(integrator,false)

and get then the error:

ERROR: LoadError: MethodError: no method matching u_modified!(::Sundials.IDA{:Dense,Nothing,Nothing}, ::Bool)
Closest candidates are:
  u_modified!(::OrdinaryDiffEq.ODEIntegrator, ::Bool) at OrdinaryDiffEq\NsugH\src\integrators\integrator_interface.jl:48
  u_modified!(::StochasticDiffEq.SDEIntegrator, ::Bool) at StochasticDiffEq\BhcRo\src\integrators\integrator_interface.jl:32
  u_modified!(::DelayDiffEq.DDEIntegrator, ::Bool) at DelayDiffEq\BCL4Q\src\integrators\interface.jl:157

I use the newest versions:

  • Sundials v4.2.5
  • DifferentialEquations v6.15.0

@ChrisRackauckas
Copy link
Member

It wasn't exported, so I just added a patch that fixes that: SciML/OrdinaryDiffEq.jl@2b3770c Should be registered in about an hour or so.

@MartinOtter
Copy link
Author

Now that NoInit() is in the docs I'll close this.

(We still need to unify Sundials with the new initialization setup, but for now it's 🤷 )

I use:

DifferentialEquations.solve(<..>, initializealg = DifferentialEquations.NoInit())   

for DifferentialEquations.IDA() solver. This works for a standard DAE. However, if an event is triggered with conditions!/affect! functions, an error is triggered when the affect! function is left:

ERROR: LoadError: Must supply differential_vars argument to DAEProblem constructor to use IDA initial value solver.

It seems NoInit() is not yet fully implemented for IDA? When supplying the differential_vars argument to DAEproblem, the event is correctly processed (but probably IDA is running its initialization code - which should be switched off).

@ChrisRackauckas
Copy link
Member

Yes it's on OrdinaryDiffEq. This is an issue for Sundials.jl then. SciML/Sundials.jl#285

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants