Hi,

I have a function that calculates the analytical solution to a pulse injection through a column (advective dispersion transport). I am trying to fit some parameters to the data.

I manage to calculate the derivative of my `pulse_injection`

function but not the derivative of the loss function. I have tried a bunch of different things without success and I was hoping someone could spot my mistake here. Thanks a lot!

This is my function:

```
function pulse_injection(
cr::Matrix{T},
x::Vector{T},
t::Vector{T},
c0,
c_in,
v,
Dl,
t_pulse
) where T
#cr = zeros(length(t), length(x))
ratio = (c0 - c_in) / 2
@inbounds for i=1:length(x), j=1:length(t)
exp_term = exp(v .* x[i] / Dl)
cr[j, i] = c_in .+ ratio .* (erfc.((x[i] .- v .* t[j])
./ (2 .* sqrt.(Dl .* t[j]))) .+ exp_term
.* erfc.((x[i] .+ v .* t[j]) ./ (2 .* sqrt.(Dl .* t[j]))))
if t[j] >= t_pulse
co = - ratio .* (erfc.((x[i] .- v .* (t[j] .- t_pulse))
./ (2 .* sqrt.(Dl .* (t[j] .- t_pulse)))) .+ exp_term
.* erfc.((x[i] .+ v .* (t[j] .- t_pulse)) ./ (2 .* sqrt.(Dl .* (t[j] .- t_pulse)))))
cr[j, i] = cr[j, i] .+ co
end
end
return nothing
end
```

This is my loss function:

```
function loss_julia(p::Vector)
phi = p[1]
alpha_l = p[2]
vp = q / phi
De = 2.01e-9 + alpha_l * vp
pulse_injection(c_model, x, t, c0, c_in, vp, De, t_pulse)
return sum((c_model[:, 1] .- c).^2)
end
```

I can run the loss function:

```
loss_julia([0.3, 1e-3])
```

and calculate the jacobian of the `pulse_injection`

function:

```
ForwardDiff.jacobian((c_model, p)-> pulse_injection(c_model, x, t, c0, c_in, q/p[1], 2.01e-9 + p[2]*q/p[1], t_pulse), c_model, [0.3, 1e-3])
```

but not from the loss function:

```
ForwardDiff.jacobian(loss_julia, [0.3, 1e-3])
```

I get this error:

```
ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_julia), Float64}, Float64, 2})```
Full code snippet:
```julia
using SpecialFunctions
using BenchmarkTools
using Plots
using Optim
using OptimizationOptimJL
using Optimization
using SciMLSensitivity
using CSV
using DataFrames
using ForwardDiff
function pulse_injection(
cr::Matrix{T},
x::Vector{T},
t::Vector{T},
c0,
c_in,
v,
Dl,
t_pulse
) where T
#cr = zeros(length(t), length(x))
ratio = (c0 - c_in) / 2
@inbounds for i=1:length(x), j=1:length(t)
exp_term = exp(v .* x[i] / Dl)
cr[j, i] = c_in .+ ratio .* (erfc.((x[i] .- v .* t[j])
./ (2 .* sqrt.(Dl .* t[j]))) .+ exp_term
.* erfc.((x[i] .+ v .* t[j]) ./ (2 .* sqrt.(Dl .* t[j]))))
if t[j] >= t_pulse
co = - ratio .* (erfc.((x[i] .- v .* (t[j] .- t_pulse))
./ (2 .* sqrt.(Dl .* (t[j] .- t_pulse)))) .+ exp_term
.* erfc.((x[i] .+ v .* (t[j] .- t_pulse)) ./ (2 .* sqrt.(Dl .* (t[j] .- t_pulse)))))
cr[j, i] = cr[j, i] .+ co
end
end
return nothing
end
using CSV
using DataFrames
#load data
data = CSV.read("data/data_a.csv",DataFrame)
global t = data[:,1]
global c = data[:,2]
global c0 = 2.0
global c_in = 0.0
global t_pulse = 3600.0
global x = [0.121]
Q0_ml = 6 #ml/hr
Q0 = Q0_ml*1e-6/3600 # m3/s
diam = 0.037 # diameter of the column [m]
area = pi*(diam/2)^2 # cross-sectional area of the column [m2]
global q = Q0/area # specific discharge [m/s]
c_model = zeros(length(t), length(x))
function loss_julia(p::Vector{T}) where T
phi = p[1]
alpha_l = p[2]
vp = q / phi
De = 2.01e-9 + alpha_l * vp
pulse_injection(c_model, x, t, c0, c_in, vp, De, t_pulse)
return sum((c_model[:, 1] .- c).^2)
end
loss_julia([0.3, 1e-3])
ForwardDiff.jacobian((c_model, p)-> pulse_injection(c_model, x, t, c0, c_in, q/p[1], 2.01e-9 + p[2]*q/p[1], t_pulse), c_model, [0.3, 1e-3])
ForwardDiff.jacobian(loss_julia, [0.3, 1e-3])
```