I am trying to solve the following problem

\begin{aligned} \nabla \cdot \nabla \phi_{1} &= \frac{\partial\left(\phi_{1}-\phi_{2}\right)}{\partial t}\\ \nabla \cdot \nabla \phi_{2} &= -\tau p_1p_2\frac{\partial\left(\phi_{1}-\phi_{2}\right)}{\partial t} \end{aligned}

in a one dimensional x\in[0, 1] domain, where p_1 and p_2 are parameters. The boundary conditions are

\begin{aligned} \phi_1(0, t)& = t \\ \nabla\phi_1(1, t) &= 0.0 \\ \phi_2(1, t) &= 0.0 \\ \nabla\phi_2(0, t) &= 0.0 \end{aligned}

I was successful in implementing the forward solver using the DAEProblem class.

using DiffEqOperators, OrdinaryDiffEq, Plots, DifferentialEquations, Zygote

nknots = 10
h = 1.0 / (nknots + 1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2

const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)

τ = 0.1
ξ = 0.5
t0 = 0.0
tlimit = 1.0 / ξ
t1 = 2.0 * tlimit
vlimit = 1.0
p = [1.0, 2.0]

function heat!(res, du, u, p, t)
rbc_1 = RobinBC((1.0, 0.0, t), (0.0, 1.0, 0.0), h, 1);
rbc_2 = RobinBC((0.0, 1.0, 0.0), (1.0, 0.0, 0.0), h, 1);
res[:, 1] = du[:, 1] - du[:, 2] - 1.0 * Δ * rbc_1 * u[:, 1]
res[:, 2] = τ * p[1] * p[2] * (du[:, 1] - du[:, 2]) + 1.0 * Δ * rbc_2 * u[:, 2]
end

u0 = zeros(nknots, 2)
diff_var = trues(nknots, 2)
du0 = zeros(nknots, 2)
prob = DAEProblem{true}(heat!, du0, u0, (t0, t1)) # , differential_vars=diff_var)
sol = solve(prob, DABDF2(), abstol=1e-3, reltol=1e-3)


Now I’d like to get gradients with respect to p_1 and p_2, with the following lines:

function sum_of_solution(u0, p)
_prob = remake(prob, u0=u0, p=p)
sum(solve(_prob, DABDF2(), rtol=1e-6, atol=1e-6, saveat=0.1, sensealg=ReverseDiffAdjoint()))
end
du01, dp1 = Zygote.gradient(sum_of_solution, u0, p)


but I am getting the error:

ERROR: LoadError: MethodError: no method matching DiffEqOperators.BoundaryPaddedVector(::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::ReverseDiff.TrackedReal{Float64, Float64, Nothing}, ::Vector{ReverseDiff.TrackedReal{Float64, Float64, ReverseDiff.TrackedArray{Float64, Float64, 2, Matrix{Float64}, Matrix{Float64}}}})
Closest candidates are:

and I think that I can only solve those problems with a DAEProblem. With that said, how can I compute the gradients? Thanks.