Error when evaluating terms on interface-coupled multifield problem

[oriolcg]

I'm getting an error when evaluating multifield components on interfaces. The error only manifests for a specific component (trial left domain against test on the right domain) and only when evaluating the residual and not the jacobian. Does anyone faced similar issue before?

Here is a reproducer of the error:

module REPRODUCER
using Gridap
model = CartesianDiscreteModel((0,1),(2,))
Ω = Interior(model)
mask = [true,false]
Ω1 = Interior(model,mask)
Ω2 = Interior(model,.!mask)
reffe = ReferenceFE(lagrangian,VectorValue{1,Float64},1)
V1 = TestFESpace(Ω1,reffe,dirichlet_tags="boundary")
V2 = TestFESpace(Ω2,reffe,dirichlet_tags="boundary")
U1 = TrialFESpace(V1,VectorValue(0.0))
U2 = TrialFESpace(V2,VectorValue(0.0))
X = MultiFieldFESpace([U1,U2])
Y = MultiFieldFESpace([V1,V2])
u1_basis = get_trial_fe_basis(U1)
u2_basis = get_trial_fe_basis(U2)
v1_basis = get_fe_basis(V1)
v2_basis = get_fe_basis(V2)
op1 = u1_basis⋅v1_basis
op2 = u2_basis⋅v1_basis
op3 = u2_basis⋅v2_basis
op4 = u1_basis⋅v2_basis
println(" op1: ", op1(Point(0.5,)))
println(" op2: ", op2(Point(0.5,)))
println(" op3: ", op3(Point(0.5,)))
println(" op4: ", op4(Point(0.5,)))
u1h = interpolate(x -> VectorValue(2x[1]),U1)
u2h = interpolate(x -> VectorValue(2x[1]),U2)
op1h = u1h⋅v1_basis
op2h = u2h⋅v1_basis
op3h = u2h⋅v2_basis
op4h = u1h⋅v2_basis
println(" op1h: ", op1h(Point(0.5,)))
println(" op2h: ", op2h(Point(0.5,)))
println(" op3h: ", op3h(Point(0.5,)))   
println(" op4h: ", op4h(Point(0.5,)))
end

And this is the error message:

op1: [0.0 0.0; 0.0 1.0]
 op2: Matrix{Float64}(undef, 2, 0)
 op3: [1.0 0.0; 0.0 0.0]
 op4: Matrix{Float64}(undef, 0, 2)
 op1h: [0.0, 1.0]
 op2h: [0.0, 0.0]
 op3h: [1.0, 0.0]
ERROR: LoadError: TypeError: in typeassert, expected Gridap.Fields.OperationField{typeof(LinearAlgebra.dot), Tuple{Gridap.Fields.LinearCombinationField{Vector{Float64}, Gridap.Fields.VoidBasis{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}, 1, Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}}, Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}}, got a value of type Gridap.Fields.OperationField{typeof(LinearAlgebra.dot), Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Fields.VoidBasis{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}, 1, Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}}, Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}}
Stacktrace:
 [1] testvalue(::Type{Gridap.Fields.OperationField{typeof(LinearAlgebra.dot), Tuple{…}}})
   @ Gridap.Fields ~/.julia/packages/Gridap/dk0DA/src/Fields/FieldsInterfaces.jl:341
 [2] testvalue(::Type{Gridap.Fields.OperationField{Gridap.Fields.OperationField{…}, Tuple{…}}})
   @ Gridap.Fields ~/.julia/packages/Gridap/dk0DA/src/Fields/FieldsInterfaces.jl:329
 [3] testitem
   @ ~/.julia/packages/Gridap/dk0DA/src/Arrays/Interface.jl:154 [inlined]
 [4] return_cache(f::Vector{Gridap.Fields.OperationField{…}}, x::Gridap.TensorValues.VectorValue{1, Float64})
   @ Gridap.Fields ~/.julia/packages/Gridap/dk0DA/src/Fields/FieldArrays.jl:5
 [5] return_cache(a::Gridap.CellData.Interpolable{…}, x::Gridap.TensorValues.VectorValue{…})
   @ Gridap.CellData ~/.julia/packages/Gridap/dk0DA/src/CellData/Interpolation.jl:36
 [6] return_cache
   @ ~/.julia/packages/Gridap/dk0DA/src/CellData/Interpolation.jl:26 [inlined]
 [7] evaluate(f::Gridap.CellData.OperationCellField{…}, x::Gridap.TensorValues.VectorValue{…})
   @ Gridap.Arrays ~/.julia/packages/Gridap/dk0DA/src/Arrays/Maps.jl:86
 [8] (::Gridap.CellData.OperationCellField{…})(x::Gridap.TensorValues.VectorValue{…})
   @ Gridap.CellData ~/.julia/packages/Gridap/dk0DA/src/CellData/CellFields.jl:247
 [9] top-level scope
   @ ~/Progs/git_repos/CMOE/Mooring.jl/reproducer.jl:36

The only difference between expected and received argument type is that one of the entries of the tree is a Vector instead of a Matrix.

[zjwegert]

Hi Oriol, I wonder if this is related to the issue I encountered in #1145?

[oriolcg]

Hi @zjwegert, yes it looks like it's the same. But I don't think it is related to AD. This reproducer is just evaluating functions on interfaces.

These are the expected and received types:

Gridap.Fields.OperationField{typeof(LinearAlgebra.dot), Tuple{Gridap.Fields.LinearCombinationField{Vector{Float64}, Gridap.Fields.VoidBasis{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}, 1, Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}}, Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}}
Gridap.Fields.OperationField{typeof(LinearAlgebra.dot), Tuple{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Fields.VoidBasis{Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}, 1, Gridap.Fields.LinearCombinationFieldVector{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}}, Gridap.Fields.LinearCombinationField{Matrix{Float64}, Gridap.Polynomials.MonomialBasis{1, Gridap.TensorValues.VectorValue{1, Float64}}}}} 

The only difference is that the first entry in the first LinearCombination is Vector{Float64} for the expected and Matrix{Float64} for the received. And cannot understand why... Any clue there?

[oriolcg]

Hi @zjwegert, this reproducer now works using this branch: https://github.com/gridap/Gridap.jl/tree/issue_%231191. Could you try to see if that also fixes the issue #1145 ?

[oriolcg]

I'm not sure if the implementation is correct and/or if it's efficient. Maybe it would be good to have a look at it @zjwegert and @JordiManyer before I plan a PR.

[oriolcg]

Looking a bit more deeper, the original error occurs because there is a LinearCombinationField that is created with values of type Vector{Float64}, while the testvalue function initializes a the LinearCombinationField with zero(eltype(V),length(fields),1), resulting in a Matrix{Float64}. Not sure why this is not encountered in other examples...

[zjwegert]

Hi @zjwegert, this reproducer now works using this branch: https://github.com/gridap/Gridap.jl/tree/issue_%231191. Could you try to see if that also fixes the issue #1145 ?

This doesn't fix #1145 because the problem in that issue is happening inside a a return_cache in FieldValues. I expect something similar to what you've done could be used to fix that.

I expect we may want to look closer at LinearCombinationField and testvalue for these two cases. Unfortunately, I won't be able to look properly at this until next week or possibly the week after.

[JordiManyer]

This fix is quite terrible performance-wise, you are instantiating the array of arrays of fields.

I expect we may want to look closer at LinearCombinationField and testvalue for these two cases.

This is what should be done. I am also currently swamped.

[oriolcg]

Another update...

The issue is that the function linear_combination creates a LinearCombinationField object from an AbstractVector in

Gridap.jl/src/Fields/FieldArrays.jl

Line 159 in 282ebae

LinearCombinationField(a,b,column)

, while the testvalue function creates a LinearCombinationField object from a matrix here:

Gridap.jl/src/Fields/FieldArrays.jl

Line 201 in 282ebae

LinearCombinationField(values,fields,1)

I now reverted the previous changes and modified the vector to a matrix in linear_combination function. Could you try again @zjwegert ?

[zjwegert]

Hi @oriolcg, I think this is on the right track! It seems to fix the typing issue in #1145, but now there is some other error

ERROR: MethodError: no method matching iterate(::VectorBlock{Vector{Int32}})
The function `iterate` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  iterate(::LibGit2.GitConfigIter)
   @ LibGit2 C:\Users\wegert\AppData\Local\Programs\Julia-1.11.6\share\julia\stdlib\v1.11\LibGit2\src\config.jl:225       
  iterate(::LibGit2.GitConfigIter, ::Any)
   @ LibGit2 C:\Users\wegert\AppData\Local\Programs\Julia-1.11.6\share\julia\stdlib\v1.11\LibGit2\src\config.jl:225       
  iterate(::Pkg.Resolve.NodePerm, Any...)
   @ Pkg C:\Users\wegert\AppData\Local\Programs\Julia-1.11.6\share\julia\stdlib\v1.11\Pkg\src\Resolve\maxsum.jl:240       
  ...

Stacktrace:
  [1] iterate
    @ .\iterators.jl:206 [inlined]
  [2] iterate
    @ .\iterators.jl:205 [inlined]
  [3] _add_entries!
    @ C:\Users\wegert\Desktop\Repos\Gridap\src\Algebra\AlgebraInterfaces.jl:163 [inlined]
  [4] add_entries!
    @ C:\Users\wegert\Desktop\Repos\Gridap\src\Algebra\AlgebraInterfaces.jl:159 [inlined]
  [5] evaluate!(cache::Nothing, k::TouchEntriesMap, A::Gridap.Algebra.CounterCSC{…}, v::Matrix{…}, i::VectorBlock{…}, j::VectorBlock{…})
    @ Gridap.Arrays C:\Users\wegert\Desktop\Repos\Gridap\src\Arrays\AlgebraMaps.jl:117
  [6] evaluate!(cache::Matrix{…}, k::TouchEntriesMap, A::Gridap.Algebra.CounterCSC{…}, v::MatrixBlock{…}, I::VectorBlock{…}, J::VectorBlock{…})
    @ Gridap.Arrays C:\Users\wegert\Desktop\Repos\Gridap\src\Arrays\AlgebraMaps.jl:148
  [7] _symbolic_loop_matrix!(A::Gridap.Algebra.CounterCSC{…}, caches::Tuple{…}, cell_rows::LazyArray{…}, cell_cols::LazyArray{…}, mat1::MatrixBlock{…}, cells::Base.OneTo{…})
    @ Gridap.FESpaces C:\Users\wegert\Desktop\Repos\Gridap\src\FESpaces\SparseMatrixAssemblers.jl:195
  [8] _symbolic_loop_matrix!(A::Gridap.Algebra.CounterCSC{…}, caches::Tuple{…}, cell_rows::LazyArray{…}, cell_cols::LazyArray{…}, mat1::MatrixBlock{…})
    @ Gridap.FESpaces C:\Users\wegert\Desktop\Repos\Gridap\src\FESpaces\SparseMatrixAssemblers.jl:190
  [9] symbolic_loop_matrix!(A::Gridap.Algebra.CounterCSC{…}, a::GenericSparseMatrixAssembler, matdata::Tuple{…})
    @ Gridap.FESpaces C:\Users\wegert\Desktop\Repos\Gridap\src\FESpaces\SparseMatrixAssemblers.jl:181
 [10] assemble_matrix(a::GenericSparseMatrixAssembler, matdata::Tuple{Vector{Any}, Vector{Any}, Vector{Any}})
    @ Gridap.FESpaces C:\Users\wegert\Desktop\Repos\Gridap\src\FESpaces\SparseMatrixAssemblers.jl:69
 [11] assemble_matrix(f::typeof(a), a::GenericSparseMatrixAssembler, U::MultiFieldFESpace{…}, V::MultiFieldFESpace{…})    
    @ Gridap.FESpaces C:\Users\wegert\Desktop\Repos\Gridap\src\FESpaces\Assemblers.jl:295
 [12] assemble_matrix(f::Function, U::MultiFieldFESpace{…}, V::MultiFieldFESpace{…})
    @ Gridap.FESpaces C:\Users\wegert\Desktop\Repos\Gridap\src\FESpaces\Assemblers.jl:364
 [13] top-level scope
    @ c:\Users\wegert\Desktop\Repos\Gridap\test.jl:14
Some type information was truncated. Use `show(err)` to see complete types.

I suspect this is something weird with doing integration on the skeleton but having one of the fields not defined.