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:
DiffEqOperators.BoundaryPaddedVector(::T, ::T, ::T2) where {T, T2<:AbstractVector{T}} at /Users/salazardetro1/.julia/packages/DiffEqOperators/FnmOt/src/boundary_padded_arrays.jl:14
I believe that my current equations might be solved using another method (maybe an ODEProblem with a mass matrix?), but my next step will be to include a nonlinear term multiplying the time derivatives as in
\begin{aligned}
\nabla \cdot \nabla \phi_{1} &= c(\phi_1, \phi_2) \frac{\partial\left(\phi_{1}-\phi_{2}\right)}{\partial t}\\
\nabla \cdot \nabla \phi_{2} &= -\tau p_1p_2c(\phi_1, \phi_2)\frac{\partial\left(\phi_{1}-\phi_{2}\right)}{\partial t}
\end{aligned}
and I think that I can only solve those problems with a DAEProblem. With that said, how can I compute the gradients? Thanks.