AutoDiff of interpolate_everywhere fails. Mismatch among FE space vector type and dual number types.

[amartinhuertas]

Hi @fverdugo,

the following MWE reproduces an error we have found when using AutoDiff (ForwardDiff) to differentiate a function that calls Gridap's interpolate_everywhere.

As shown in the below, the error can be workarounded using interpolate_everywhere! and vectors of DoFs of the appropriate eltype (i.e.. dual number types).

The cause of the error is a mismatch among the FE space vector type (Float64) and the result of evaluating the function to be interpolated with dual numbers.

Do you see a neat approach to solve this issue? I guess that an ad-hoc differentiation rule for interpolate_everywhere that calls interpolate_everywhere! under the hood to compute the derivative might be one possible approach. Anything else that come intos your mind? E.g. creating a FE space with dual numbers as vector type with e.g., gradient(::FESpace) ?

Thanks!

using Gridap
using FillArrays
using Test

function dP_dp_error(params)
  function g(params)
    f(x)=(params[1]+params[2]*params[2]*x[1]+(params[3]^3)*x[2])
    VectorValue(get_free_dof_values(
     (Gridap.FESpaces.interpolate_everywhere(f,P))))
  end
  gradient(g)(VectorValue(params))
end

function dP_dp_workaround(params)
  function g(params)
    f(x)=(params[1]+params[2]*params[2]*x[1]+(params[3]^3)*x[2])
    free_values = Vector{eltype(params)}(undef,num_free_dofs(P))
    dirichlet_values = Vector{eltype(params)}(undef,num_dirichlet_dofs(P))
    VectorValue(get_free_dof_values(
     (Gridap.FESpaces.interpolate_everywhere!(f, free_values,dirichlet_values,P))))
  end
  gradient(g)(VectorValue(params))
end

domain = (0,1,0,1)
ncells = (5,5)
model = CartesianDiscreteModel(domain,ncells)

# FE space to represent parameter space
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
P  = TestFESpace(model,reffe; conformity=:H1)

params=[1.0,2.0,3.0]
dP_dp_params=dP_dp_workaround(params)
dP_dp_params=dP_dp_error(params)

[amartinhuertas]

cc @wei3li @santiagobadia

[fverdugo]

@amartinhuertas IMO the workaround already looks better than anything else I can think of.

[amartinhuertas]

Ok. So we should find a way to automatically trigger the workaround whenever autodiff tries to differentiate interpolate_everywhere. Correct?

[fverdugo]

Ok. So we should find a way to automatically trigger the workaround whenever autodiff tries to differentiate interpolate_everywhere. Correct?

Or perhaps use the workaround directly. Since I guess it is not easy to write a rule for interpolate_everywhere. The variables you differentiate wrt are not passed explicitly to interpolate_everywhere, but captured in the given clousure...

If you find a way to define this rule, nice, but perhaps it is not worth the effort since the workaround is just a couple on lines longer than the targeted code.

[amartinhuertas]

Or perhaps use the workaround directly. Since I guess it is not easy to write a rule for interpolate_everywhere. The variables you differentiate wrt are not passed explicitly to interpolate_everywhere, but captured in the given clousure...

Ok. I understand what you mean. We also though about this. We can define a new type, say, ParameterizedFunction (instead of a closure) to make the parameters explicit.

[amartinhuertas]

Sorry ... ParameterizedFunction, NOT ParameterizedFEFunction, solved in the previous post.