AutoDiff on sub-domains

[ConnorMallon]

There are some cases when using Gridap and Gridap's Autodiff machinery together does not work. Here we have a problem where a FEFunction only exists in certain parts of the domain. It means that if we do something like get_cell_dof_values we get back an array for which some of the entries are empty. We have problems then using Autodiff in the multifield case where one field is empty and the other isnt at a given point. The following MWE gives the target functionality:

module AutodiffSubDomainTests

using Test
using Gridap
import Gridap: ∇
using LinearAlgebra: tr, ⋅
using Gridap
import Gridap: ∇
using Test
using Gridap.FESpaces
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.MultiField

# Background model
n = 10
mesh = (n,n)
domain = (0,1,0,1) .- 1
order = 1
model = CartesianDiscreteModel(domain, mesh)
Ω = Triangulation(model)

# Extract a portion of the background mesh
R = 0.7
function is_in(coords)
  n = length(coords)
  x = (1/n)*sum(coords)
  d = x[1]^2 + x[2]^2 - R^2
  d < 0
end
cell_to_coods = get_cell_coordinates(Ω)
cell_to_is_in = collect1d(lazy_map(is_in,cell_to_coods))
cell_to_is_out = collect1d(lazy_map(!,cell_to_is_in))

Ω_in = Triangulation(model,cell_to_is_in)
Ω_out = Triangulation(model,cell_to_is_out)

dΩ_in = Measure(Ω_in,2)
dΩ_out = Measure(Ω_out,2)

reffe = ReferenceFE(lagrangian,Float64,order)
V_in  = TestFESpace(Ω_in, reffe)
V_out = TestFESpace(Ω_out,reffe)

V = MultiFieldFESpace([V_in,V_out])

uh = FEFunction(V,ones(num_free_dofs(V)))
j((u1,u2)) = ∫(u1)dΩ_in + ∫(u2)dΩ_out
js=∇(j)(uh)
assemble_vector(js,V)

end # module

[fverdugo]

Hi @ConnorMallon I can reproduce this in the release_0.17 branch with the following error:

ERROR: LoadError: AssertionError: A check failed
Stacktrace:
 [1] macro expansion
   @ ~/Code/jl/Gridap.jl/src/Helpers/Macros.jl:60 [inlined]
 [2] get_triangulation(f::Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain})
   @ Gridap.MultiField ~/Code/jl/Gridap.jl/src/MultiField/MultiFieldCellFields.jl:34
 [3] get_triangulation
   @ ~/Code/jl/Gridap.jl/src/MultiField/MultiFieldFEFunctions.jl:30 [inlined]
 [4] _gradient(f::Function, uh::Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain}}, fuh::Gridap.CellData.DomainContribution)
   @ Gridap.FESpaces ~/Code/jl/Gridap.jl/src/FESpaces/FEAutodiff.jl:22
 [5] gradient(f::typeof(Main.Issue661.j), uh::Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain}})
   @ Gridap.FESpaces ~/Code/jl/Gridap.jl/src/FESpaces/FEAutodiff.jl:5
 [6] (::Gridap.Fields.var"#_gradient#83"{typeof(Main.Issue661.j)})(x::Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain}})
   @ Gridap.Fields ~/Code/jl/Gridap.jl/src/Fields/AutoDiff.jl:20
 [7] top-level scope
   @ ~/Code/jl/Gridap.jl/issue_661.jl:49
 [8] include(fname::String)
   @ Base.MainInclude ./client.jl:444
 [9] top-level scope
   @ REPL[2]:1
in expression starting at /home/fverdugo/Code/jl/Gridap.jl/issue_661.jl:1

automatic differentiation for multi-field is not fully implemented now.

[fverdugo]

Unfortunately this is not a dummy bug that can be solved easily. It is a whole functionality still not implemented in the library.

To implement this one needs to understand:

• a bit the internals of ForwardDiff
• how AD wrt FEFunctions is implemented in Gridap
• how FE spaces are generated in a portion of the domain (as of the new version in release_0.17)
• how mulit-field is implemented in Gridap, specifically the case in which the single-fields are not defined on the same domain

On top of this, we first want to solve issue #584 since the solution of this issue can affect on how AD is implemented for multi-field

[ConnorMallon]

Hey @amartinhuertas @fverdugo, seeing as multifield autodiff is fixed it would be nice to also have this functionality. The first problem I am having is in the change_argument function:

function FESpaces._change_argument(
  op::typeof(gradient),f,trian,uh::MultiFieldFEFunction)
  U = get_fe_space(uh)
  function g(cell_u)
    single_fields = GenericCellField[]
    nfields = length(U.spaces)
    cell_dofs_field_offsets=_get_cell_dofs_field_offsets(uh)
    for i in 1:nfields
      view_range=cell_dofs_field_offsets[i]:cell_dofs_field_offsets[i+1]-1
      cell_values_field = lazy_map(a->view(a,view_range),cell_u)
      cf = CellField(U.spaces[i],cell_values_field)
      push!(single_fields,cf)
    end
    xh = MultiFieldCellField(single_fields)
    cell_grad = f(xh)
    get_contribution(cell_grad,trian)
  end
  g
end

The context is when we call the function in a specific sub-domain (trian). The _get_cell_dofs_field_offsets(uh) can be specialised for the given trian _get_cell_dofs_field_offsets(uh,trian) easily but I am not sure how to specialise the building of the cellfield cf = CellField(U.spaces[i],cell_values_field) for the trian. At the moment the cell_basis from the space U.spaces[i] is not defined on the same trian as the cell_values_field. How can we make these the same ? I guess there should be a change_domain function for dealing with this case ?

Please see a5fd715 for what I have done in my forked repo and the target functionality test to reproduce the error.

[oriolcg]

Hi @ConnorMallon , @amartinhuertas , @fverdugo, I'm getting the following error when trying to use autodiff for a problem with mixed-dimensional PDEs (a variable defined in a 2D domain and another defined on one of the boundaries):

ERROR: LoadError: AssertionError: 

This method does not make sense for multi-field
since each field can be defined on a different triangulation.
Pass a triangulation in the second argument to get the DOF values
on top of the corresponding cells.

Stacktrace:
 [1] macro expansion at C:\Users\ocolomesgene\.julia\packages\Gridap\mK4Q1\src\Helpers\Macros.jl:60 [inlined]
 [2] get_cell_dof_values(::Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain}}) at C:\Users\ocolomesgene\.julia\packages\Gridap\mK4Q1\src\MultiField\MultiFieldFEFunctions.jl:44
 [3] _jacobian(::Function, ::Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain}}, ::Gridap.CellData.DomainContribution) at C:\Users\ocolomesgene\.julia\packages\Gridap\mK4Q1\src\MultiField\MultiFieldFEAutodiff.jl:46
 [4] jacobian(::GridapODEs.TransientFETools.var"#res_0#82"{Float64,GridapODEs.TransientFETools.TransientCellField{Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain}}},Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain},typeof(Main.FreeSurfacePotentialFlowTests.res)}, ::Gridap.MultiField.MultiFieldFEFunction{Gridap.MultiField.MultiFieldCellField{Gridap.CellData.ReferenceDomain}}) at C:\Users\ocolomesgene\.julia\packages\Gridap\mK4Q1\src\FESpaces\FEAutodiff.jl:49

Is this error related to this issue? If not, is it easier to solve?

[ConnorMallon]

Hey @oriolcg, that problem can be solved by adding the triangulation that you are looping over at that point as an argument to the get_cell_dof_values function (see LOC of fork ) but it will not solve the whole problem since there are also other parts of the _jacobain,_gradient etc. that should also depend on the triangulations which do not. I think this issue would need to be solved before using autodiff on any problem which involves fields defined on different triangulation so i guess it is the same issue.

[amartinhuertas]

@oriolcg @ConnorMallon I will take a look at it whenever I have some time. Thanks for reporting

[amartinhuertas]

@oriolcg do you have a mwe?

[oriolcg]

I'll prepare one reproducer and post it here

[ConnorMallon]

hey @amartinhuertas, have you had a chance to take a look at this?

[amartinhuertas]

I'll prepare one reproducer and post it here

Hi @oriolcg ... could you please help us and prepare one reproducer of this issue? I want to work with both @ConnorMallon reproducers and yours, so that I can be sure that I am finding a solution that considers both cases (if they are diferent)

[oriolcg]

Hi @amartinhuertas , yes I forgot about this issue. Here I attach a reproducer:

module tmp

using Gridap

domain = (0.0, 1.0, 0.0, 1.0)
partition = (3,3)
model_Ω = CartesianDiscreteModel(domain,partition)
labels = get_face_labeling(model_Ω)
add_tag_from_tags!(labels,"surface",[3,4,6])
Ω = Interior(model_Ω)
Γ = Boundary(model_Ω,tags="surface")
order = 1
degree = 2*order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)
reffe = ReferenceFE(lagrangian,Float64,order)
V_Ω = TestFESpace(Ω,reffe,conformity=:H1)
V_Γ = TestFESpace(Γ,reffe,conformity=:H1)
U_Ω = TrialFESpace(V_Ω)
U_Γ = TrialFESpace(V_Γ)
Y = MultiFieldFESpace([V_Ω,V_Γ])
X = MultiFieldFESpace([U_Ω,U_Γ])
res((ϕ,η),(w,v)) = ∫(∇(ϕ)⋅∇(w) )dΩ + ∫( η*v )dΓ
op = FEOperator(res,X,Y)
(ϕₕ,ηₕ) = solve(op)
println(∑(∫(ϕₕ*ϕₕ)dΩ))
println(∑(∫(ηₕ*ηₕ)dΓ))

end

[amartinhuertas]

Hi @amartinhuertas , yes I forgot about this issue. Here I attach a reproducer:

Thanks!

[oriolcg]

Hi @ConnorMallon @amartinhuertas , what's the status of this issue. I am still facing problems using autodiff in multifield problems defined in different domains. Do you have any work-around for this?

[ConnorMallon]

Hey @oriolcg, I ended up not needing this functionality so I don't have a work-around. I am still happy however to help this get implemented.