added support of complex numbers for mesh adaptivity in MIRK method

[JulianUlrich]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Complex states seemed not to be compatible with mesh adaptivity. I just took the magnitude of the values in cache instead of their complex representation. Maybe not the best way, but it works...

[github-actions]

Benchmark Results

	master	c87aef1...	master/c87aef1493fabe...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2()	6.65 ± 0.16 ms	6.45 ± 0.18 ms	1.03
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3()	2.44 ± 0.12 ms	2.44 ± 0.11 ms	1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4()	0.863 ± 0.052 ms	0.856 ± 0.055 ms	1.01
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5()	0.905 ± 0.072 ms	0.904 ± 0.066 ms	1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6()	1.02 ± 0.1 ms	1.02 ± 0.11 ms	1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false)	1.98 ± 0.57 ms	1.96 ± 0.58 ms	1.01
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true)	3.58 ± 0.49 ms	3.62 ± 0.83 ms	0.989
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.0561 ± 0.0083 s	0.0565 ± 0.0075 s	0.992
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.0677 ± 0.0093 s	0.0719 ± 0.014 s	0.942
Simple Pendulum/IIP/Shooting(Tsit5())	0.256 ± 0.077 ms	0.256 ± 0.077 ms	1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2()	0.0412 ± 0.0029 s	0.0395 ± 0.003 s	1.04
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3()	11.9 ± 0.82 ms	11.6 ± 0.66 ms	1.02
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4()	3.41 ± 0.21 ms	3.27 ± 0.24 ms	1.04
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5()	3.49 ± 0.27 ms	3.35 ± 0.26 ms	1.04
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6()	3.58 ± 0.25 ms	3.49 ± 0.29 ms	1.03
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false)	3.49 ± 2.6 ms	3.59 ± 2.6 ms	0.972
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true)	6.15 ± 4.6 ms	6.36 ± 4.7 ms	0.967
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false)	0.103 ± 0.0051 s	0.108 ± 0.0052 s	0.957
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true)	0.126 ± 0.0093 s	0.134 ± 0.013 s	0.939
Simple Pendulum/OOP/Shooting(Tsit5())	0.612 ± 0.044 ms	0.612 ± 0.046 ms	1
time_to_load	7.09 ± 0.054 s	7.02 ± 0.015 s	1.01

Benchmark Plots

A plot of the benchmark results have been uploaded as an artifact to the workflow run for this PR.
Go to "Actions"->"Benchmark a pull request"->[the most recent run]->"Artifacts" (at the bottom).

[ErikQQY]

@JulianUlrich Thank you for your contribution! Under what circumstances do you want to use this?

[JulianUlrich]

I got an error from adaptivity when having complex states in the model. This can be useful in electrical engineering. F.e. for evaluating the impedance of a transmission line model. Mostly, only the boundary values are needed as a result, so it makes sense to start with a coarse mesh.

[ErikQQY]

Yeah, I can understand using complex states when we have impedances in the circuit, but I suspect the current implementation can support such models, can you please share an MWE or some simple models that can illustrate the basic idea, there may be more things to do to support these cases.

[JulianUlrich]

Here's a simple code for testing:

using BoundaryValueDiffEq
using Plots

xspan = (0.0,1);

function simpleTLM!(du,u,p,t)
    Φ  = u[1];
    dΦ = u[2];
    j = 0.1 * Φ + 1im*Φ;
    du[1] = dΦ;
    du[2] = 10*(1-1*t)*j;
end

function bc2a!(resid_a, u_a, p)
    resid_a[1] = u_a[2] + Complex{Float64}(1)
end

function bc2b!(resid_b, u_b, p)
    resid_b[1] = u_b[2] 
end

bvp2 = TwoPointBVProblem(simpleTLM!, (bc2a!, bc2b!), [Complex{Float64}(0),Complex{Float64}(1)], xspan;
    bcresid_prototype = (zeros(ComplexF64,1), zeros(ComplexF64,1)))
@time sol2 = solve(bvp2, MIRK4(), dt = 1, reltol=1e-9)

x = mapreduce(permutedims, vcat, sol2.u)

pl = plot(sol2.t, abs.(x))
display(pl)

This could be a Transmission line with space-dependent conductivity...
When using the code without the changes to adaptivity.jl, i get the following error because of the complex states:

ERROR: MethodError: no method matching isless(::Int64, ::ComplexF64)
The function `isless` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  isless(::Missing, ::Any)
   @ Base missing.jl:87
  isless(::Any, ::Missing)
   @ Base missing.jl:88
  isless(::Integer, ::ForwardDiff.Dual{Ty}) where Ty
   @ ForwardDiff user\.julia\packages\ForwardDiff\UBbGT\src\dual.jl:149
  ...

Stacktrace:
  [1] <(x::Int64, y::ComplexF64)
    @ Base .\operators.jl:353
  [2] mesh_selector!(cache::BoundaryValueDiffEqMIRK.MIRKCache{…})
    @ BoundaryValueDiffEqMIRK user\.julia\packages\BoundaryValueDiffEqMIRK\wGmgf\src\adaptivity.jl:51

[ErikQQY]

Yeah, the main issue here is the adaptivity routines, could you please add this as a test case and also add the changes to the FIRK methods adaptivity as well?