SplitODE is not differentiable?

[rveltz]

Hi!

I would like to obtain the derivative of the flow at fixed time. However, it errors with

ERROR: MethodError: no method matching exponential!(::Base.ReshapedArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#17#18", Float64}, Float64, 1}, 2, SubArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#17#18", Float64}, Float64, 1}, 1, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#17#18", Float64}, Float64, 1}}, Tuple{UnitRange{Int64}}, true}, Tuple{}}, ::ExponentialUtilities.ExpMethodHigham2005Base)
Closest candidates are:
  exponential!(::StridedMatrix{T}, ::ExponentialUtilities.ExpMethodHigham2005Base) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at ~/.julia/packages/ExponentialUtilities/e1lsW/src/exp_baseexp.jl:24
  exponential!(::StridedMatrix{T}, ::ExponentialUtilities.ExpMethodHigham2005Base, ::Any) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at ~/.julia/packages/ExponentialUtilities/e1lsW/src/exp_baseexp.jl:24
  exponential!(::Any) at ~/.julia/packages/ExponentialUtilities/e1lsW/src/exp.jl:15
  ...
Stacktrace:
  [1] phiv_dense!(w::Base.ReshapedArray{ForwardDiff.Dual{ForwardDiff.Tag{var"#17#18", Float64}, Float6

MWE:

using Revise
	using DiffEqOperators, ForwardDiff, DifferentialEquations, SparseArrays
	using LinearAlgebra, Plots, Parameters, Setfield

function Laplacian2D(Nx, Ny, lx, ly, bc = :Dirichlet)
	hx = 2lx/Nx
	hy = 2ly/Ny
	D2x = CenteredDifference(2, 2, hx, Nx)
	D2y = CenteredDifference(2, 2, hy, Ny)
	if bc == :Neumann
		Qx = Neumann0BC(hx)
		Qy = Neumann0BC(hy)
	elseif  bc == :Dirichlet
		Qx = Dirichlet0BC(typeof(hx))
		Qy = Dirichlet0BC(typeof(hy))
	end
	D2xsp = sparse(D2x * Qx)[1]
	D2ysp = sparse(D2y * Qy)[1]

	A = kron(sparse(I, Ny, Ny), D2xsp) + kron(D2ysp, sparse(I, Nx, Nx))
	return A, D2x
end

@inbounds function NL!(f, u, p, t = 0.)
	@unpack r, μ, ν, c3, c5 = p
	n = div(length(u), 2)
	u1 = @view u[1:n]
	u2 = @view u[n+1:2n]

	f1 = @view f[1:n]
	f2 = @view f[n+1:2n]

	for i=1:n
		ua = u1[i]^2 + u2[i]^2
		f1[i] = r * u1[i] - ν * u2[i] - ua * (c3 * u1[i] - μ * u2[i]) - c5 * ua^2 * u1[i]
		f2[i] = r * u2[i] + ν * u1[i] - ua * (c3 * u2[i] + μ * u1[i]) - c5 * ua^2 * u2[i]
	end

	return f
end

####################################################################################################
Nx = 41*1
	Ny = 21*1
	n = Nx*Ny
	lx = pi
	ly = pi/2

	Δ = Laplacian2D(Nx, Ny, lx, ly)[1]
	par_cgl = (r = 1.2, μ = 0.1, ν = 1.0, c3 = -1.0, c5 = 1.0, Δ = blockdiag(Δ, Δ))
	sol0 = 0.1rand(2Nx, Ny)
	sol0_f = vec(sol0)

f1 = DiffEqArrayOperator(par_cgl.Δ)
f2 = NL!
prob_sp = SplitODEProblem(f1, f2, sol0_f, (0.0, 1.0), @set par_cgl.r = 1.2; reltol = 1e-8, dt = 0.1)
sol = @time solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1)

function flow(x0, prob0)
	prob = remake(prob0, u0 = x0)
	sol = solve(prob, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1)
	return sol[end]
end

flow(sol0_f, prob_sp)

ForwardDiff.derivative(t->flow(sol0_f .+ t .* sol0_f, prob_sp), zero(eltype(sol0_f)))

[rveltz]

Hi,

is it an easy fix or a mistake on my side?

[ChrisRackauckas]

Extending that type signature beyond StridedMatrix may be all that's required.