Refactor of autodiff for MultiField #1070
Draft
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1070 +/- ##
==========================================
+ Coverage 89.05% 89.24% +0.19%
==========================================
Files 197 197
Lines 23905 23767 -138
==========================================
- Hits 21288 21212 -76
+ Misses 2617 2555 -62 ☔ View full report in Codecov by Sentry. |
|
All the GridapTopOpt state map tests pass for this branch. I've also benchmarked the changes and compared to master. Here are the results for 6e649c5:
Scriptusing Gridap, Gridap.MultiField
using BenchmarkTools
function main(n)
println("Running: n = $n")
model = CartesianDiscreteModel((0,1,0,1),(n,n))
order = 1
Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)
V1 = FESpace(Ω,ReferenceFE(lagrangian,Float64,order);dirichlet_tags="boundary")
V2 = FESpace(Ω,ReferenceFE(lagrangian,VectorValue{2,Float64},order);dirichlet_tags="boundary")
U1 = TrialFESpace(V1)
U2 = TrialFESpace(V2)
UB = MultiFieldFESpace([U1,U2])
VB = MultiFieldFESpace([V1,V2])
uh = FEFunction(UB,rand(num_free_dofs(UB)))
res((u1,u2),(v1,v2)) = ∫(u1*∇(u1)⋅∇(v1) + (u2 ⋅ u2) * ∇(u2)⊙∇(v2))dΩ
## Gradient
println("Gradient:")
_res_grad(uh) = res(uh,uh)
g = gradient(_res_grad,uh)
@btime gradient($_res_grad,$uh);
@btime Gridap.FESpaces.assemble_vector($g,$VB);
## Jacobian
println("Jacobian:")
dv = get_fe_basis(VB)
_res_jac(uh) = res(uh,dv)
j = jacobian(_res_jac,uh)
@btime jacobian($_res_jac,$uh);
@btime Gridap.FESpaces.assemble_matrix($j,$UB,$VB);
## Hessian
println("Hessian:")
h = hessian(_res_grad,uh);
@btime hessian($_res_grad,$uh);
@btime Gridap.FESpaces.assemble_matrix($h,$UB,$VB);
end
main(64);
main(128);
main(256); |
|
Note that this currently breaks for the following using Gridap
order = 1
degree = 2order
model = CartesianDiscreteModel((0,1,0,1),(10,10))
V_φ = TestFESpace(model,ReferenceFE(lagrangian,Float64,1))
φf(x) = sqrt((x[1]-0.5)^2 + (x[2]-0.5)^2) - 0.3
φh = interpolate(φf,V_φ)
V = TestFESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},order,space=:P))
Q = TestFESpace(model,ReferenceFE(lagrangian,Float64,order,space=:P))
VQ = MultiFieldFESpace([V,Q])
Λ = SkeletonTriangulation(model)
dΛ = Measure(Λ,degree)
n_Λ = get_normal_vector(Λ)
dΩ = Measure(get_triangulation(model),degree)
res((u,p),(v,q),φ) = ∫(u⋅v + φ)dΩ + ∫(jump(n_Λ ⋅ ∇(u)) ⋅ jump(n_Λ ⋅ ∇(v)))dΛ
# res((u,p),(v,q),φ) = ∫(u⋅v + φ)dΩ + ∫(jump(n_Λ ⋅ ∇(u)) ⋅ jump(n_Λ ⋅ ∇(v)) + 0mean(φ))dΛ # workaround
∇(φ->res(zero(VQ),zero(VQ),φ),φh)with ERROR: MethodError: no method matching +(::Int64, ::Tuple{Int64, Int64})
The function `+` exists, but no method is defined for this combination of argument types.
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any...)
@ Base operators.jl:596
+(::Real, ::Complex{Bool})
@ Base complex.jl:322
+(::Integer, ::BlockArrays.Block)
@ BlockArrays ~/.julia/packages/BlockArrays/X84bj/src/blockindices.jl:73
...
Stacktrace:
[1] extract_gradient_block!
@ ~/.julia/dev/Gridap/src/Fields/ArrayBlocks.jl:1583 [inlined]
[2] evaluate!
@ ~/.julia/dev/Gridap/src/Fields/ArrayBlocks.jl:1536 [inlined]
[3] evaluate
@ ~/.julia/dev/Gridap/src/Arrays/Maps.jl:87 [inlined]
[4] return_value
@ ~/.julia/dev/Gridap/src/Arrays/Maps.jl:64 [inlined]
[5] return_type
@ ~/.julia/dev/Gridap/src/Arrays/Maps.jl:62 [inlined]
[6] lazy_map(::Gridap.Arrays.AutoDiffMap, ::Gridap.Arrays.LazyArray{…}, ::Gridap.Arrays.LazyArray{…})
@ Gridap.Arrays ~/.julia/dev/Gridap/src/Arrays/LazyArrays.jl:57
[7] autodiff_array_gradient(a::Gridap.FESpaces.var"#g#66"{…}, i_to_x::Gridap.Arrays.LazyArray{…}, j_to_i::Gridap.Arrays.LazyArray{…})
@ Gridap.Arrays ~/.julia/dev/Gridap/src/Arrays/Autodiff.jl:31
[8] _gradient(f::Function, uh::Gridap.FESpaces.SingleFieldFEFunction{…}, fuh::Gridap.CellData.DomainContribution)
@ Gridap.FESpaces ~/.julia/dev/Gridap/src/FESpaces/FEAutodiff.jl:23
[9] gradient(f::var"#15#16", uh::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{ReferenceDomain}})
@ Gridap.FESpaces ~/.julia/dev/Gridap/src/FESpaces/FEAutodiff.jl:5
[10] top-level scope
@ ~/.julia/dev/GridapTopOpt/scripts/Embedded/Bugs/multifield_MWE.jl:23
Some type information was truncated. Use `show(err)` to see complete types.The current workaround is included in the script |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
The current implementation of autodiff for MultiField spaces is a mess... It relies on a lot of specific code, and has severe performance issues. I believe that with some though we could use the same code as for SingleField.
This is motivated by the fact that I do not want to replicate the current implementation for GridapDistributed. I believe it will also solve many issues we currently have for autodiff+multifield.
To-Dos
To-Think
Skeletons are wrong when the number of dofs is different in the left and right cells. The issue is that we are dualizing both left and right cells at once, but witout knowing which cells those are (we just assume they are the same). I think I'll have to go back to what Kishore was doing, i.e doing left then right. Hopefully without any densifying maps.
The other solution would be to directly port the cell values to the destination triangulation, and evaluate the bilinear forms with cellfields that are already there. The issue is that this will likely break the change_domain apis in an unpredictable way.