Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Computing gradient wrt Dirichlet data #969

Open
arnold-pdev opened this issue Dec 27, 2023 · 0 comments
Open

Computing gradient wrt Dirichlet data #969

arnold-pdev opened this issue Dec 27, 2023 · 0 comments

Comments

@arnold-pdev
Copy link

arnold-pdev commented Dec 27, 2023

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/[email protected]/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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant