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)
The error I receive is
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}})
@ Main ~/Library/CloudStorage/GoogleDrive-andrewarnold@math.arizona.edu/Shared drives/Atmosphere 2/regression/LapFunctionBox.jl:54
[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]
[13] vector_mode_gradient(f::var"#377#378", x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:89
[14] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}}, ::Val{true})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:0
[15] gradient(f::Function, x::Vector{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#377#378", Float64}, Float64, 10}}})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:17
[16] gradient(f::Function, x::Vector{Float64})
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/gradient.jl:17
[17] top-level scope
@ REPL[111]:1`
Clearly at some point the Dirichlet data is being defined as Float64, but I hope you might know where and how to resolve this issue. I’m using Optim.jl to solve this optimization problem, but if you know a way to solve this problem within the Gridap.jl ecosystem, I’m all ears.