# Computing sensitivity to Dirichlet data in Gridap.jl

Hi all, happy holidays-

I’m trying to use Gridap to invert the Laplace problem to determine Dirichlet data on the boundaries given some point measurements on the interior. The idea is to minimize the squared residual between simulated and measured data. However, I am having trouble using ForwardDiff (finite differencing works, but I really want AD to speed up the optimization process).

Here is my program:

``````# Finite element method
using Gridap
using GridapGmsh, Gridap.Io
model = DiscreteModelFromFile("model03-1e4.json")
order = 2
reffe = ReferenceFE(lagrangian,Float64,order)
# define the test space
V0 = TestFESpace(model,reffe;conformity=:H1,
dirichlet_tags=[
"westFace",
"eastFace",
"southFace",
"northFace",
"westBase",
"eastBase",
"southBase",
"northBase",
"upFace",
"bottom"]
)
degree = 2
Ω = Triangulation(model)
dΩ = Measure(Ω,degree) # approximation of Lebesgue measure

function laplacesolve(α)
"""
laplacesolve solves the Laplace equation with Dirichlet data given by the vector α.
"""
# add boundary data (here: fixed data for boundary regions 1:10)
g_wf(x) = α[1]  # west face
g_ef(x) = α[2]  # east face
g_sf(x) = α[3]  # south face
g_nf(x) = α[4]  # north face
g_wb(x) = α[5]  # west base
g_eb(x) = α[6]  # east base
g_sb(x) = α[7]  # south base
g_nb(x) = α[8]  # north base
g_uf(x) = α[9]  # up face
g_bo(x) = α[10] # bottom

# define the trial space
Ug = TrialFESpace(V0,[g_wf,g_ef,g_sf,g_nf,g_wb,g_eb,g_sb,g_nb,g_uf,g_bo])

# define the weak problem
f(x) = 0.0 # ⇒ Laplace
a(u,v) =∫( ∇(u)⋅∇(v) )*dΩ # Dirichlet terms are built into the TrialFESpace, not here
b(v) = ∫( f*v )*dΩ #+ ∫( v*h )*dΓ #for a neumann term
# define the problem through the abstract type FEOperator
op = AffineFEOperator(a,b,Ug,V0)
ls = LUSolver()
solver = LinearFESolver(ls)
uh = solve(solver,op) # solve

# extract node values at pts
return uh
end

function tempout(bcs::T, pts) where T
uh = laplacesolve(bcs)
# extract node values at pts
return uh.(Point.(pts))
end
``````

I’m trying to differentiate a function (bcs)->f(bcs,data,pts), where f is defined as
`f(bcs,data,pts) = norm(tempout(bcs,pts) .- data, 2)^2/length(data)`

``````ERROR: MethodError: no method matching
Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10})

Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:207
(::Type{T})(::T) where T<:Number
@ Core boot.jl:792
(::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
@ Base char.jl:50`
``````

### Stacktrace

``````Stacktrace:
[1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10})
@ Base ./number.jl:7
[2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}, i1::Int64)
@ Base ./array.jl:969
[3] setindex!
@ ./multidimensional.jl:670 [inlined]
[4] _free_and_dirichlet_values_fill!(free_vals::Vector{Float64}, dirichlet_vals::Vector{Float64}, cache_vals::Tuple{Tuple{Nothing, Tuple{Gridap.Arrays.CachedVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}, Tuple{Gridap.Arrays.CachedVector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}, Nothing, Tuple{Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}}}}, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Nothing, Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, TensorValue{3, 3, Float64, 9}}}, Tuple{Tuple{Nothing, Nothing, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{VectorValue{3, Float64}, Vector{VectorValue{3, Float64}}}, Tuple{Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{VectorValue{3, Float64}}}}, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, VectorValue{3, Float64}}}}}, Gridap.Arrays.IndexItemPair{Int64, Gridap.Fields.AffineMap{3, 3, Float64, 9}}}}}, Gridap.Arrays.IndexItemPair{Int64, Gridap.Fields.OperationField{Gridap.Fields.GenericField{var"#g_wf#340"{Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}}, Tuple{Gridap.Fields.AffineMap{3, 3, Float64, 9}}}}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}}, cache_dofs::Gridap.Arrays.CachedVector{Int32, Vector{Int32}}, cell_vals::Gridap.Arrays.LazyArray{Gridap.Arrays.CompressedArray{Gridap.ReferenceFEs.LagrangianDofBasis{VectorValue{3, Float64}, Int64}, 1, Vector{Gridap.ReferenceFEs.LagrangianDofBasis{VectorValue{3, Float64}, Int64}}, Vector{Int8}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{typeof(∘)}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.OperationField{Gridap.Fields.GenericField{var"#g_wf#340"{Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}}, Tuple{Gridap.Fields.AffineMap{3, 3, Float64, 9}}}, 1, Tuple{FillArrays.Fill{Gridap.Fields.GenericField{var"#g_wf#340"{Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(Gridap.Fields.affine_map), 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, VectorValue{3, Float64}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{Float64}, 1, Vector{Vector{Float64}}, Vector{Int8}}}}}}}}}}, cell_dofs::Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}, cells::Vector{Int32})
@ Gridap.FESpaces ~/.julia/packages/Gridap/I60jU/src/FESpaces/UnconstrainedFESpaces.jl:135
[5] gather_dirichlet_values!(dirichlet_vals::Vector{Float64}, f::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Nothing}, cell_vals::Gridap.Arrays.LazyArray{Gridap.Arrays.CompressedArray{Gridap.ReferenceFEs.LagrangianDofBasis{VectorValue{3, Float64}, Int64}, 1, Vector{Gridap.ReferenceFEs.LagrangianDofBasis{VectorValue{3, Float64}, Int64}}, Vector{Int8}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{typeof(∘)}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.OperationField{Gridap.Fields.GenericField{var"#g_wf#340"{Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}}, Tuple{Gridap.Fields.AffineMap{3, 3, Float64, 9}}}, 1, Tuple{FillArrays.Fill{Gridap.Fields.GenericField{var"#g_wf#340"{Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{typeof(Gridap.Fields.affine_map), 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.AffineMap{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, TensorValue{3, 3, Float64, 9}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{VectorValue{3, Float64}}, 1, Vector{Vector{VectorValue{3, Float64}}}, Vector{Int8}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.LinearCombinationMap{Colon}, 1, Tuple{Base.OneTo{Int64}}}, VectorValue{3, Float64}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Reindex{Vector{VectorValue{3, Float64}}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{VectorValue{3, Float64}}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.CompressedArray{Vector{Float64}, 1, Vector{Vector{Float64}}, Vector{Int8}}}}}}}}}})
@ Gridap.FESpaces ~/.julia/packages/Gridap/I60jU/src/FESpaces/UnconstrainedFESpaces.jl:106
[6] compute_dirichlet_values_for_tags!(dirichlet_values::Vector{Float64}, dirichlet_values_scratch::Vector{Float64}, f::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Nothing}, tag_to_object::Vector{Function})
@ Gridap.FESpaces ~/.julia/packages/Gridap/I60jU/src/FESpaces/SingleFieldFESpaces.jl:250
[7] compute_dirichlet_values_for_tags(f::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Nothing}, tag_to_object::Vector{Function})
@ Gridap.FESpaces ~/.julia/packages/Gridap/I60jU/src/FESpaces/SingleFieldFESpaces.jl:237
[8] TrialFESpace(space::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Nothing}, objects::Vector{Function})
@ Gridap.FESpaces ~/.julia/packages/Gridap/I60jU/src/FESpaces/TrialFESpaces.jl:23
[9] laplacesolve(α::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}})
[10] tempout(bcs::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}, pts::Vector{VectorValue{3, Float64}})
@ Main ./REPL[110]:2
[11] (::var"#377#378")(bcs::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}})
@ Main ./REPL[111]:1
[12] vector_mode_dual_eval!
@ ~/.julia/packages/ForwardDiff/PcZ48/src/apiutils.jl:24 [inlined]