MultiFieldFESpace with complex numbers

[simonsticko]

Hi.

I tried to use MultiFieldFESpace based on two complex spaces, but I am getting a type conversion error. Is this a bug? Or have I misunderstood how this should be done? (Gridap v0.17.22)

Minimal working example:

using Gridap

# Using float works fine. Complex gives an error.
# T = Float64
T = ComplexF64

domain = (0, 1, 0, 1)
partition = (4, 4)
model = CartesianDiscreteModel(domain, partition)

order = 1
reffe = ReferenceFE(lagrangian, Float64, order)

V1 = TestFESpace(model, reffe; conformity=:H1, vector_type=Vector{T})
V2 = TestFESpace(model, reffe; conformity=:H1, vector_type=Vector{T})

U1 = TrialFESpace(V1)
U2 = TrialFESpace(V2)

Y = MultiFieldFESpace([V1, V2])
X = MultiFieldFESpace([U1, U2])

degree = 2 * order
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

# Project constant 1 into both spaces.
a((u1, u2), (v1, v2)) = ∫(v1 * u1)dΩ + ∫(v2 * u2)dΩ
l((v1, v2)) = ∫(v1 + v2)dΩ

op = AffineFEOperator(a, l, X, Y)
uh1, uh2 = solve(op)

Error message:

ERROR: LoadError: MethodError: Cannot `convert` an object of type 
  Gridap.Fields.ArrayBlock{Array{Float64,1},1} to an object of type 
  Gridap.Fields.ArrayBlock{Array{ComplexF64,1},1}

Closest candidates are:
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:84

Stacktrace:
  [1] cvt1
    @ ./essentials.jl:463 [inlined]
  [2] ntuple
    @ ./ntuple.jl:49 [inlined]
  [3] convert(::Type{Tuple{Gridap.Fields.MatrixBlock{Matrix{Float64}}, Gridap.Fields.VectorBlock{Vector{ComplexF64}}}}, x::Tuple{Gridap.Fields.MatrixBlock{Matrix{Float64}}, Gridap.Fields.VectorBlock{Vector{Float64}}})
    @ Base ./essentials.jl:465
  [4] setproperty!
    @ ./Base.jl:40 [inlined]
  [5] getindex!(cache::Tuple{Tuple{Nothing, Tuple{Gridap.Fields.VectorBlock{Gridap.Arrays.CachedVector{ComplexF64, Vector{ComplexF64}}}, Tuple{Gridap.Fields.VectorBlock{Vector{ComplexF64}}, Vector{Nothing}}}, Tuple{Nothing, Tuple{Tuple{Nothing, Gridap.Fields.VectorBlock{Vector{ComplexF64}}, Tuple{Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{ComplexF64, Vector{ComplexF64}}, Tuple{Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{ComplexF64}}}, Tuple{Tuple{Nothing, Gridap.Arrays.CachedVector{ComplexF64, Vector{ComplexF64}}, Tuple{Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}}, Gridap.Arrays.IndexItemPair{Int64, Vector{ComplexF64}}}}}, Gridap.Arrays.IndexItemPair{Int64, Gridap.Fields.VectorBlock{Vector{ComplexF64}}}}, Tuple{Tuple{Nothing, Nothing, Tuple{Nothing, Nothing}}, Gridap.Arrays.IndexItemPair{Int64, Bool}}}}, Gridap.Arrays.IndexItemPair{Int64, Tuple{Gridap.Fields.MatrixBlock{Matrix{Float64}}, Gridap.Fields.VectorBlock{Vector{ComplexF64}}}}}, a::Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.CellData.AttachDirichletMap, 1, Tuple{Base.OneTo{Int64}}}, Tuple{Gridap.Fields.MatrixBlock{Matrix{Float64}}, Gridap.Fields.VectorBlock{Vector{ComplexF64}}}, 1, Tuple{FillArrays.Fill{Tuple{Gridap.Fields.MatrixBlock{Matrix{Float64}}, Gridap.Fields.VectorBlock{Vector{Float64}}}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Fields.BlockMap{1}, 1, Tuple{Base.OneTo{Int64}}}, Gridap.Fields.VectorBlock{Vector{ComplexF64}}, 1, Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Gridap.Arrays.PosNegReindex{Gridap.Arrays.SubVector{ComplexF64, Vector{ComplexF64}}, Vector{ComplexF64}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{ComplexF64}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Broadcasting{Gridap.Arrays.PosNegReindex{Gridap.Arrays.SubVector{ComplexF64, Vector{ComplexF64}}, Vector{ComplexF64}}}, 1, Tuple{Base.OneTo{Int64}}}, Vector{ComplexF64}, 1, Tuple{Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}}}}}, Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.MultiField.var"#34#36", 1, Tuple{Base.OneTo{Int64}}}, Bool, 1, Tuple{Vector{Bool}, Vector{Bool}}}}}, i::Int64)
    @ Gridap.Arrays ~/.julia/packages/Gridap/JQVke/src/Arrays/LazyArrays.jl:214
  [6] numeric_loop_matrix_and_vector!(A::Gridap.Algebra.InserterCSC{ComplexF64, Int64}, b::Vector{ComplexF64}, a::Gridap.FESpaces.GenericSparseMatrixAssembler, data::Tuple{Tuple{Vector{Any}, Vector{Any}, Vector{Any}}, Tuple{Vector{Any}, Vector{Any}, Vector{Any}}, Tuple{Vector{Any}, Vector{Any}}})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JQVke/src/FESpaces/SparseMatrixAssemblers.jl:311
  [7] assemble_matrix_and_vector(a::Gridap.FESpaces.GenericSparseMatrixAssembler, data::Tuple{Tuple{Vector{Any}, Vector{Any}, Vector{Any}}, Tuple{Vector{Any}, Vector{Any}, Vector{Any}}, Tuple{Vector{Any}, Vector{Any}}})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JQVke/src/FESpaces/SparseMatrixAssemblers.jl:101
  [8] AffineFEOperator(weakform::Gridap.FESpaces.var"#60#61"{typeof(a), typeof(l)}, trial::MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{ComplexF64}}, test::MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{ComplexF64}}, assem::Gridap.FESpaces.GenericSparseMatrixAssembler)
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JQVke/src/FESpaces/AffineFEOperators.jl:38
  [9] AffineFEOperator(::Function, ::MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{ComplexF64}}, ::Vararg{MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{ComplexF64}}})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JQVke/src/FESpaces/AffineFEOperators.jl:47
 [10] AffineFEOperator(::typeof(a), ::typeof(l), ::MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{ComplexF64}}, ::Vararg{MultiFieldFESpace{Gridap.MultiField.ConsecutiveMultiFieldStyle, Gridap.FESpaces.UnConstrained, Vector{ComplexF64}}})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/JQVke/src/FESpaces/AffineFEOperators.jl:51
 [11] top-level scope
    @ ~/workspace/gridap-bug/test.jl:31
in expression starting at /home/simon/workspace/gridap-bug/test.jl:31

[JordiManyer]

Hey @simonsticko, I've managed to reproduce the error.

For some reason, the code is not allocating the elemental contributions correctly (creating float-based arrays instead of complex-based arrays).

A complete fix might be a little involved, but for now you can do this:

using Gridap

# Using float works fine. Complex gives an error.
# T = Float64
T = ComplexF64

domain = (0, 1, 0, 1)
partition = (4, 4)
model = CartesianDiscreteModel(domain, partition)

order = 1
reffe = ReferenceFE(lagrangian, Float64, order)

V1 = TestFESpace(model, reffe; conformity=:H1, vector_type=Vector{T})
V2 = TestFESpace(model, reffe; conformity=:H1, vector_type=Vector{T})

U1 = TrialFESpace(V1)
U2 = TrialFESpace(V2)

Y = MultiFieldFESpace([V1, V2])
X = MultiFieldFESpace([U1, U2])

degree = 2 * order
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

ONE = 1.0+0.0im

# Project constant 1 into both spaces.
a((u1, u2), (v1, v2)) = ∫(ONE*(v1 * u1 + v2 * u2))dΩ
l((v1, v2)) = ∫(ONE*(v1 + v2))dΩ

op = AffineFEOperator(a, l, X, Y)
uh1, uh2 = solve(op)

It's basically your code, but I force the type of the integrand to be of complex type by multiplying everywhere by a complex-valued one.

Actually: your weakform has nothing complex about it, it only involves real-valued functions. When you try an actual complex-valued problem, complex numbers will appear naturally and you will no longer have to do this trick (at least that would be my guess).

Let me know if it works for you as well.

[simonsticko]

Let me know if it works for you as well.

Yes. That works. Thank you for the work-around, and for the quick reply. 👍

When you try an actual complex-valued problem, complex numbers will appear naturally and you will no longer have to do this trick (at least that would be my guess).

I don't think this is necessarily true. In for example, Helmholz equation the coefficients in the weak form are real, unless there is damping added.