AD for multifield and skeletons

[JordiManyer]

Sister PR for this.

[zjwegert]

This fails when calling AD for multifield when the spaces are created from the triangulation instead of the model:

function main_mf_Ω(distribute,parts)
  ranks = distribute(LinearIndices((prod(parts),)))

  model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(4,4))

  k = 2
  reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},k)
  reffe_p = ReferenceFE(lagrangian,Float64,k-1;space=:P)

  Ω = Triangulation(model)
  dΩ = Measure(Ω,2*(k+1))

  u(x) = VectorValue(x[2],-x[1])
  V = TestFESpace(Ω,reffe_u,dirichlet_tags="boundary")            # Taking Ω instead of model
  U = TrialFESpace(V,u)
  Q = TestFESpace(Ω,reffe_p;conformity=:L2,constraint=:zeromean)  # Taking Ω instead of model

  X = MultiFieldFESpace([U,Q])
  Y = MultiFieldFESpace([V,Q])

  ν = 1.0
  f = VectorValue(0.0,0.0)

  conv(u,∇u) = (∇u')⋅u
  dconv(du,∇du,u,∇u) = conv(u,∇du)+conv(du,∇u)
  c(u,v) = ∫(v⊙(conv∘(u,∇(u))))dΩ
  dc(u,du,dv) = ∫(dv⊙(dconv∘(du,∇(du),u,∇(u))))dΩ

  biform((du,dp),(dv,dq)) = ∫(ν*∇(dv)⊙∇(du) - (∇⋅dv)*dp - (∇⋅du)*dq)dΩ
  liform((dv,dq)) = ∫(dv⋅f)dΩ

  r((u,p),(v,q)) = biform((u,p),(v,q)) + c(u,v) - liform((v,q))
  j((u,p),(du,dp),(dv,dq)) = biform((du,dp),(dv,dq)) + dc(u,du,dv)

  op = FEOperator(r,j,X,Y)
  op_AD = FEOperator(r,X,Y)

  xh = interpolate([VectorValue(1.0,1.0),1.0],X)
  uh, ph = xh
  A = jacobian(op,xh)
  A_AD = jacobian(op_AD,xh)
  @test reduce(&,map(≈,partition(A),partition(A_AD)))

  g((v,q)) = ∫(0.5*v⋅v + 0.5*q*q)dΩ
  dg((v,q)) = ∫(uh⋅v + ph*q)dΩ
  b = assemble_vector(dg,X)
  b_AD = assemble_vector(gradient(g,xh),X)
  @test b ≈ b_AD
end

A (likely suboptimal) fix for now is something like

function test_triangulation(Ω1,Ω2)
  @assert typeof(Ω1.grid) == typeof(Ω2.grid)
  t = map(fieldnames(typeof(Ω1.grid))) do field
    getfield(Ω1.grid,field) == getfield(Ω2.grid,field)
  end
  all(t)
  a = Ω1.model === Ω2.model
  b = Ω1.tface_to_mface == Ω2.tface_to_mface
  a && b && all(t)
end

function CellData.get_triangulation(f::MultiFieldCellField)
  s1 = first(f.single_fields)
  trian = get_triangulation(s1)
  # @check all(map(i->trian===get_triangulation(i),f.single_fields))
  @check all(map(i->test_triangulation(trian,get_triangulation(i)),f.single_fields))
  trian
end

Note that CellData.get_triangulation(a::DistributedMultiFieldCellField) fails as well because the triangulations are not identical objects - a similar fix to the above can fix this as well.

I suspect the reason this occurs is because of the first line in

function FESpaces.FESpace( _trian::DistributedTriangulation,reffe;split_own_and_ghost=false,constraint=nothing,kwargs...)
  trian = add_ghost_cells(_trian) # <-----
  spaces = map(local_views(trian)) do t
    FESpace(t,reffe;kwargs...)
  end
  gids = generate_gids(trian,spaces)
  vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost)
  space = DistributedSingleFieldFESpace(spaces,gids,trian,vector_type)
  return _add_distributed_constraint(space,reffe,constraint)
end