Shape derivative

[ConnorMallon]

I am trying to make GridapEmbedded.jl compatible with ForwardDiff.jl to the point in which we could autodiff an integral ∫( f )dΩ with respect to parameters describing the domain e.g. level set values:

d( ∫( f )dΩ ) / d(level_set_values) where the domain we wish to integrate over depends on the level set values

To use ForwardDiff's autodiff machinery, we must allow the function that maps level set values to the integral ∫( a(u,v) )dΩ to also accept as an input the type ForwardDiff.Dual{} <: Real. Type restrictions in Gridap and GridapEmbedded mean that julia tries to promote the Dual{} type to a Float64 giving an error.

I have tried to remove some of the concrete type restrictions in Gridap / GridapEmbedded but have only managed as to get so far as being able to compute d( ∫( f )dΩ ) / d(level_set_values) where f is a function of the domain e.g. a constant. See commits: https://github.com/ConnorMallon/Gridap.jl/commit/16df3a0835804d5e41ce4d1d69c87bc2fff88145 https://github.com/ConnorMallon/GridapEmbedded.jl/commit/54560e22aa9d4c83fc9822676490dbeee40a2874

We would however want to be able to compute the derivative when f is composed of cellfields/ fefunctions. I am stuck here as I can not find where in the code the types are restricted. A minimum example of what we would want working:

module ShapeDiff

using ForwardDiff
using GridapEmbedded
using GridapEmbedded.LevelSetCutters
using Gridap
using Gridap.ReferenceFEs
using Gridap.Arrays

domain = (0,1,0,1); cells=(10,10)
bgmodel = CartesianDiscreteModel(domain,cells)

point_to_coords = collect1d(get_node_coordinates(bgmodel))
geo1 = disk(0.3,x0=Point(0.5,0.5))
geo1_d = discretize(geo1,bgmodel) 
lvl_set = geo1_d.tree.data[1]

geo = DiscreteGeometry(lvl_set,point_to_coords,name="")
cutgeo = cut(bgmodel,geo)
model = DiscreteModel(cutgeo)
order=1
Vstd = FESpace(model,FiniteElements(PhysicalDomain(),model,lagrangian,Float64,order))
u0 = interpolate(1,Vstd)

function residual(lvl_set)
  geo = DiscreteGeometry(lvl_set,point_to_coords,name="")
  cutgeo = cut(bgmodel,geo)
  Ω = Triangulation(cutgeo)
  dΩ = Measure(Ω,2)
  #sum(∫(1)dΩ)
  sum(∫(u0)dΩ)
end

residual(lvl_set)

drdl(lvl_set) =  ForwardDiff.gradient(residual,lvl_set)
@show drdl(lvl_set)

end 

[fverdugo]

Hi @ConnorMallon

as you noticed, one needs to propagate dual numbers through the code. This implies, to be less strict in some places and replace e.g. Float64 by something more general. Note however, that replacing e.g. Float64 by e.g. Real is too naive since it would introduce a lot of type in-stabilities in the code. The solution is to add more type parameters, when required. E.g.,

Replace

struct CutTriangulation{Dc,A}
  sub_trian::A
  done_ls_to_cell_to_inoutcut::Vector{Vector{Int8}}
  pending_ls_to_point_to_value::Vector{Vector{Float64}}
  minicell::MiniCell
  table::LookupTable{Dc}
end

by

struct CutTriangulation{Dc,A,T}
  sub_trian::A
  done_ls_to_cell_to_inoutcut::Vector{Vector{Int8}}
  pending_ls_to_point_to_value::Vector{Vector{T}}
  minicell::MiniCell
  table::LookupTable{Dc}
end

instead of

struct CutTriangulation{Dc,A}
  sub_trian::A
  done_ls_to_cell_to_inoutcut::Vector{Vector{Int8}}
  pending_ls_to_point_to_value::AbstractVector
  minicell::MiniCell
  table::LookupTable{Dc}
end

as I can see in your code.

[ConnorMallon]

Hey @fverdugo,

I have updated the generalisations using type parameters as you say: 87e3dc05779e9c7f58acfb8943bb45b39e3ab7b7. I still have the problem that I cannot trace back to where some of the objects get fixed to Float64's in Gridap. I have only found a few as you can see that have been updated in the commits but there are still some remaining which I cannot find which are preventing me from performing the differentiation. Do you have any idea where they could be or have any tips for helping me to find them?

[ConnorMallon]

I think i found it - I will submit a PR soon.