Complex woes with implicit solvers (out-of-place only)

[antoine-levitt]

using OrdinaryDiffEq, DiffEqBase

fun(dgamma, gamma,p,t) = (dgamma .= imag(gamma))
gamma0 = [1.0+im 0;0 1]
# gamma0 = eye(2,2)
prob = ODEProblem(fun,gamma0,(0.,1.))
sol = solve(prob,Trapezoid(autodiff=false))

Works for Rosenbrock23, not for Trapezoid or other solvers in the sdirk family. It fails on the else branch of

  if alg.tol != nothing
    tol = alg.tol
  else
    tol = min(0.03,first(reltol)^(0.5))
  end

saying that it can't compare real and complex. Specifying a real reltol gives another error that the constructor does not match the structure. I think some type needs a real() somewhere.

[ChrisRackauckas]

Yeah, thanks for the issue. Before we had complex support in the differentiation tools this wasn't an issue, but now we really need to go back through and test complex numbers everywhere.

[ChrisRackauckas]

Yeah, that tolerance might have to be real-valued.

[MSeeker1340]

I looked into this a bit and the problem might be more complex than I originally thought.

For most integrators that support complex-valued systems, adding a real conversion inside the init method for integrators should be enough of a fix:

# In src/solve.jl:113
  if typeof(alg) <: FunctionMap
    abstol_internal = real.(zero(u))
  elseif abstol == nothing
    if uBottomEltypeNoUnits == uBottomEltype || !(typeof(u) <: ArrayPartition)
      abstol_internal = real(uBottomEltype(uBottomEltype(1)*1//10^6))
    else
      abstol_internal = real.(ones(u).*1//10^6)
    end
  else
    abstol_internal = real.(abstol)
  end

  if typeof(alg) <: FunctionMap
    reltol_internal = real.(zero(first(u)/t))
  elseif reltol == nothing
    reltol_internal = real(uBottomEltypeNoUnits(1//10^3))
  else
    reltol_internal = real.(reltol)
  end

However for the sdirk as well as kencarp/kverno family of solvers, the situation is more complex because they store custom tolerance values and related data inside their caches. Take the implicit Euler for example:

# in src/caches/sdirk_caches.jl:1
mutable struct ImplicitEulerCache{uType,rateType,uNoUnitsType,J,UF,JC,uEltypeNoUnits,F} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::J
  W::J
  uf::UF
  jac_config::JC
  linsolve::F
  ηold::uEltypeNoUnits
  κ::uEltypeNoUnits
  tol::uEltypeNoUnits
  newton_iters::Int
end

# in src/caches/sdirk_caches.jl:42
  
# in src/caches/sdirk_caches.jl:52
    tol = min(0.03,first(reltol)^(0.5))

ηold, κ and tol should all be real numbers, so their type specifications are wrong if uEltypeNoUnits is complex. To fix this, we can include an additional type parameter (say uToltype) in ImplicitEulerCache, change all real parameters to this type, and fix statements like ηold = one(uEltypeNoUnits) in the cache and perform_step functions. Repeat this for all methods in the two families and we should be fine.

I think a more systematic approach should nevertheless be considered because the above process is very error-prone. Even methods outside these families are probably not 100% safe, because a comparison statement may lurk somewhere that will cause complex systems to fail from trying to compare two complex values.

Also, the fact that these two families of methods store the tolerance information inside the caches instead of integrator.opts.abstol and integrator.opts.reltol to be a bit inconsistent. It may be worthwhile to modify init so that the tolerance options can be set more intelligently.

[antoine-levitt]

I think a more systematic approach should nevertheless be considered because the above process is very error-prone. Even methods outside these families are probably not 100% safe, because a comparison statement may lurk somewhere that will cause complex systems to fail from trying to compare two complex values.

I have a vested interest in making sure that the julia ecosystem in general is complex numbers-friendly, so if somebody does the hard work of fixing the solvers, I can write tests to ensure that doesn't happen. (although it certainly looks as if sdirk_cache could use a good refactoring)

[ChrisRackauckas]

although it certainly looks as if sdirk_cache could use a good refactoring

Yes, it certainly can.

Also, the fact that these two families of methods store the tolerance information inside the caches instead of integrator.opts.abstol and integrator.opts.reltol to be a bit inconsistent. It may be worthwhile to modify init so that the tolerance options can be set more intelligently.

init only has the common interface arguments.

http://docs.juliadiffeq.org/latest/basics/common_solver_opts.html

Otherwise if every algorithm could add new arguments to it, there would be a bunch of minor arguments that are only used by one obscure non-diagonal SDE solver (that example would be true). So instead algorithms can have algorithm specific arguments and (error-checks) as part of their algorithm constructor and this allows them to document it more fully like you see with the Sundials algorithms.

http://docs.juliadiffeq.org/latest/solvers/ode_solve.html#Sundials.jl-1

These tolerances for the SDIRK solvers are algorithm specific and only apply to this family. While there can definitely be some refactoring to make this area more concise (I'd go even further and say it's badly needed) and to make this arg-handling less error-prone, making this argument available everywhere doesn't truly make sense because it's very specific to this form of quasi-Newton iteration for a family of methods in this package's implementation.

[ChrisRackauckas]

I have a vested interest in making sure that the julia ecosystem in general is complex numbers-friendly, so if somebody does the hard work of fixing the solvers, I can write tests to ensure that doesn't happen.

Would you mind submitting some tests and adding @test_broken for now? That would be helpful.

The biggest hurdle to complex support in stiff solvers has been the lack of differencing tools. Now that DiffEqDiffTools.jl is the backend and can handle that, we need to go back through the stiff solvers and find out which ones need to be changed to support it.

Additionally, we can get Sundials.jl supporting complex numbers by, in the common interface, handling the conversion to and from complex numbers to higher dimensional arrays (just via reinterpret). That project shouldn't be all that bad if someone wants to give it a try.

[MSeeker1340]

Since I have already gone through the code for the caches for the SDIRK family, let me do an ad hoc fix using the uToltype approach so that the solvers can handle complex systems and leave the refactoring for later.

[ChrisRackauckas]

The inplace solvers were fixed up. Out of place needs more work and a larger refactor. Thanks @MSeeker1340

[ChrisRackauckas]

SciML/OrdinaryDiffEq.jl@f612cfd