Complex Matrix Support

[WenjieYao]

These lines unconditionally creates a Float64 matrix which fails if my weak form returns complex numbers. What is the good way to fix this?

[stevengj]

It would be nice to have a way early on to specify the desired number field relatively early on and have this information propagate through subsequent steps.

(For example, conceivably one might want not only complex matrices, but also dimensionful matrices using Unitful.jl, or quaternion matrices, double-double matrices, etcetera. Though obviously complex matrices will be the most common need beyond Float64.)

[stevengj]

(We need this for a Helmholtz equation with perfectly matched layer (PML) boundary layers, which lead to complex matrix entries.)

[fverdugo]

Gridap is quite agnostic to the eltype, but not completely as you say. This is an old issue that we are aware of.

For reference, could you include here some snipped of the weak form you are willing to compute (or a simplified version that has the same problem)?

[fverdugo]

@stevengj @WenjieYao

BTW, the fe operators and assemblers have optional arguments to specify matrix and vector types. E.g.

https://github.com/gridap/Gridap.jl/blob/4c94645c2cb44e60d7efd48de01f8826616957bf/src/FESpaces/AffineFEOperators.jl#L51

you could do some think like

mtype = SparseMatrixCSC{Complex{Float64},Int,Int}
op = AffineFEOperator(mtype,U,V,terms...)

(even though I have never tried myself with complex numbers)

[WenjieYao]

For a scalar magnetic Helmholtz equation, the weak form looks like this: $a(u,v)=\frac{1}{\epsilon}\nabla u\cdot\nabla v-k_0^2\mu u\cdot v$, it will return complex value if $\epsilon$ is complex (for lossy materials) or a PML is implemented (The PML is equivalent to a material of a complex $\epsilon$).

[fverdugo]

@WenjieYao

For a scalar magnetic Helmholtz equation, the weak form looks like this: $a(u,v)=\frac{1}{\epsilon}\nabla u\cdot\nabla v-k_0^2\mu u\cdot v$, it will return complex value if $\epsilon$ is complex (for lossy materials) or a PML is implemented (The PML is equivalent to a material of a complex $\epsilon$).

OK, then try to use

mtype = SparseMatrixCSC{Complex{Float64},Int}
op = AffineFEOperator(mtype,U,V,terms...)

and see if it works

[fverdugo]

I also assume that you will need to extract the real part of the solution once you have solved the system, right?

[WenjieYao]

I can successfully create a AffineFEOperator op with mytype, but then the problem is with the solve(op), it will still give an InexactError: Float64(156.3643713837483 + 686.7830771274798im). Although I can still solve the equation with some low-level functions that can return the matrix A,b, it would be better if the solve function also supports complex numbers.

For reference, here's the stack trace

 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:847 [inlined]
 [4] _unsafe_copyto!(::Array{Float64,1}, ::Int64, ::Array{Complex{Float64},1}, ::Int64, ::Int64) at ./array.jl:257
 [5] unsafe_copyto! at ./array.jl:311 [inlined]
 [6] _copyto_impl! at ./array.jl:335 [inlined]
 [7] copyto! at ./array.jl:321 [inlined]
 [8] copyto! at ./array.jl:347 [inlined]
 [9] copy_entries! at /Users/jayyao/.julia/packages/Gridap/p256b/src/Algebra/AlgebraInterfaces.jl:57 [inlined]
 [10] solve! at /Users/jayyao/.julia/packages/Gridap/p256b/src/Algebra/LinearSolvers.jl:255 [inlined]
 [11] solve! at /Users/jayyao/.julia/packages/Gridap/p256b/src/Algebra/LinearSolvers.jl:193 [inlined]
 [12] solve!(::Array{Float64,1}, ::LUSolver, ::Gridap.Algebra.AffineOperator{SparseMatrixCSC{Complex{Float64},Int64},Array{Float64,1}}) at /Users/jayyao/.julia/packages/Gridap/p256b/src/Algebra/NonlinearSolvers.jl:22
 [13] solve!(::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{true,()}}, ::LinearFESolver, ::AffineFEOperator, ::Nothing) at /Users/jayyao/.julia/packages/Gridap/p256b/src/FESpaces/FESolvers.jl:103
 [14] solve!(::Gridap.FESpaces.SingleFieldFEFunction{Gridap.CellData.GenericCellField{true,()}}, ::LinearFESolver, ::AffineFEOperator) at /Users/jayyao/.julia/packages/Gridap/p256b/src/FESpaces/FESolvers.jl:15
 [15] solve(::LinearFESolver, ::AffineFEOperator) at /Users/jayyao/.julia/packages/Gridap/p256b/src/FESpaces/FESolvers.jl:37
 [16] solve(::AffineFEOperator) at /Users/jayyao/.julia/packages/Gridap/p256b/src/FESpaces/FESolvers.jl:43
 [17] top-level scope at In[4]:71
 [18] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091```

[WenjieYao]

For the solution, I need both real and imaginary part

[fverdugo]

It seems that you need a complex valued fe space. In any case, in order to be able to help further i need the precise mathematical definition of the weak form you are trying to solve, including the definition of the discrete test and trial fe spaces.

[stevengj]

So we should use valuetype=ComplexF64 when defining the test and trial functions?

(Even though we have a complex-valued quadratic form a(u,v)=∫[ε⁻¹(x)∇u⋅∇v - k²μ(x)uv], because the spatially varying coefficients ε and μ are complex numbers, and hence complex-valued test and trial functions u,v, the actual FE basis functions don't need to be complex-valued… only the coefficients need to be complex.)

[stevengj]

Setting the valuetype=ComplexF64 gives a MethodError in the TestFESpace constructor.

[stevengj]

@WenjieYao will create a self-contained "toy version" of the problem and post it here.

[WenjieYao]

Here's a toy version of the problem:

using Gridap
using SparseArrays
domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

order = 1
V0 = TestFESpace(
  reffe=:Lagrangian, order=order, valuetype=ComplexF64,
  conformity=:H1, model=model, dirichlet_tags="boundary")

u(x) = x[1] + x[2]
f(x) = 0
U = TrialFESpace(V0,u)

trian = Triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)

fϵ(x) = 1.0+0.1im
a(u,v) = fϵ*∇(v)⊙∇(u)
b(v) = v*f

t_Ω = AffineFETerm(a,b,trian,quad)
matrix_type = SparseMatrixCSC{ComplexF64,Int}
vector_type = Vector{ComplexF64}
op = AffineFEOperator(matrix_type,vector_type,U,V0,t_Ω)

#uh = solve(op)

[fverdugo]

Setting the valuetype=ComplexF64 gives a MethodError in the TestFESpace constructor.

Yes, this is known. Fixing this should be relatively easy. I Hope I can find some time after an important deadline to investigate this.

[stevengj]

Thanks! (In the meantime, we have a workaround by calling lower-level functions.)

[oriolcg]

Has this issue been solved?

[fverdugo]

Yes!

You can have real-valued shape functions and complex dof values, to generate a complex-valued interpolations.

https://github.com/gridap/Gridap.jl/blob/d1236879510c6703221a036926a004782ad75546/test/FESpacesTests/SingleFieldFESpacesTests.jl#L43

[zlin-opt]

Hi, is this issue solved? I am trying to generate a complex PML matrix for linear elastodynamic PDE where you need to combine a complex tensor S with the symmetric gradient like this:

S = TensorValue{3,3,ComplexF64}(1+1im,0,0,0,1+1im,0,0,0,1+1im) εs(u) = (S⋅outer(∇,u) + outer(u,∇)⋅S)/2 σ(εs) = λ tr(εs) one(εs) + 2 μ εs

But it is giving me InexactError: Float64(0.0 + 1.4555555555555542e10im)

I have included a minimal code below; thanks for any help.

using Gridap domain = (0,10,0,10,0,1) partition = (20,20,5) model = CartesianDiscreteModel(domain, partition) order=1 degree=2*order

Ω = Triangulation(model) dΩ = Measure(Ω,degree)

reffe = ReferenceFE(lagrangian,VectorValue{3,ComplexF64},order) V = TestFESpace(model, reffe; conformity=:H1) U = V

S = TensorValue{3,3,ComplexF64}(1+1im,0,0,0,1+1im,0,0,0,1+1im)

λ=5e10 μ=3e10 σ(εs) = λ tr(εs) one(εs) + 2 μ εs εs(u) = (S⋅outer(∇,u) + outer(u,∇)⋅S)/2

Λ(u,v) = εs(v)⊙(σ∘εs(u)) A_mat = assemble_matrix(U,V) do u, v ∫(Λ(u,v))dΩ end

[amartinhuertas]

Try to use ComplexF64 type all the way through, even for those scalars which are mathematically real-valued scalars as the Lame constants.

[zlin-opt]

Try to use ComplexF64 type all the way through, even for those scalars which are mathematically real-valued scalars as the Lame constants.

I changed: λ=5e10 + 0im μ=3e10 + 0im

but still the same error.

[amartinhuertas]

Do it for all scalars, not only the Lame constants

[zlin-opt]

Do it for all scalars, not only the Lame constants

I even changed S = TensorValue{3,3,ComplexF64}(1+1im,0+0im,0+0im,0+0im,1+1im,0+0im,0+0im,0+0im,1+1im) Is there any other scalar I missed? I can't seem to find anything else than λ and μ σ εs, Λ are functions which supposedly should be type flexible

[amartinhuertas]

There are two 2s

[zlin-opt]

There are two 2s

Changed them: σ(εs) = λ tr(εs) one(εs) + (2+0im) μ εs εs(u) = (S⋅outer(∇,u) + outer(u,∇)⋅S)/(2+0im)

but still the same error.

[amartinhuertas]

Ok.

At this stage, what I would do is to start simplifying the bilinear form progressively to see if there is a particular term that is responsible for the error.

For example, define lambda as:

Λ(u,v) = εs(v)⊙εs(u)

and simplify:

σ(εs) = μ * εs

In any case, do you have the full stack of the error?

[zlin-opt]

Ok.

At this stage, what I would do is to start simplifying the bilinear form progressively to see if there is a particular term that is responsible for the error.

For example, define lambda as:

Λ(u,v) = εs(v)⊙εs(u)

and simplify:

σ(εs) = μ * εs

In any case, do you have the full stack of the error?

I simplified everything to just:

εs(u) = S⋅outer(∇,u) Λ(u,v) = εs(v)⊙εs(u)

Still the same error. So I guess the problem is with the outer, ⋅ or ∇ or ⊙ Here is the error stack:

InexactError: Float64(0.0 + 0.36666666666666625im)

Stacktrace: [1] Real @ .\complex.jl:44 [inlined] [2] convert @ .\number.jl:7 [inlined] [3] setindex! @ .\array.jl:966 [inlined] [4] add_entry!(#unused#::typeof(+), a::Gridap.Algebra.InserterCSC{Float64, Int64}, v::ComplexF64, i::Int32, j::Int32) @ Gridap.Algebra C:\Users\zinli.julia\packages\Gridap\Qj36b\src\Algebra\SparseMatrixCSC.jl:132 [5] _add_entries! @ C:\Users\zinli.julia\packages\Gridap\Qj36b\src\Algebra\AlgebraInterfaces.jl:149 [inlined] [6] add_entries! @ C:\Users\zinli.julia\packages\Gridap\Qj36b\src\Algebra\AlgebraInterfaces.jl:127 [inlined] [7] evaluate!(cache::Nothing, k::Gridap.Arrays.AddEntriesMap{typeof(+)}, A::Gridap.Algebra.InserterCSC{Float64, Int64}, v::Matrix{ComplexF64}, i::Vector{Int32}, j::Vector{Int32}) @ Gridap.Arrays C:\Users\zinli.julia\packages\Gridap\Qj36b\src\Arrays\AlgebraMaps.jl:47 [8] _numeric_loop_matrix!(mat::Gridap.Algebra.InserterCSC{Float64, Int64}, caches::Tuple{Nothing, Nothing, Gridap.Arrays.CachedVector{Int32, Vector{Int32}}, Gridap.Arrays.CachedVector{Int32, Vector{Int32}}}, cell_vals::FillArrays.Fill{Matrix{ComplexF64}, 1, Tuple{Base.OneTo{Int64}}}, cell_rows::Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}, cell_cols::Gridap.Arrays.Table{Int32, Vector{Int32}, Vector{Int32}}) @ Gridap.FESpaces C:\Users\zinli.julia\packages\Gridap\Qj36b\src\FESpaces\SparseMatrixAssemblers.jl:227 [9] numeric_loop_matrix!(A::Gridap.Algebra.InserterCSC{Float64, Int64}, a::Gridap.FESpaces.GenericSparseMatrixAssembler, matdata::Tuple{Vector{Any}, Vector{Any}, Vector{Any}}) @ Gridap.FESpaces C:\Users\zinli.julia\packages\Gridap\Qj36b\src\FESpaces\SparseMatrixAssemblers.jl:214 [10] assemble_matrix(a::Gridap.FESpaces.GenericSparseMatrixAssembler, matdata::Tuple{Vector{Any}, Vector{Any}, Vector{Any}}) @ Gridap.FESpaces C:\Users\zinli.julia\packages\Gridap\Qj36b\src\FESpaces\SparseMatrixAssemblers.jl:71 [11] assemble_matrix(f::var"#1#2", a::Gridap.FESpaces.GenericSparseMatrixAssembler, U::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Gridap.FESpaces.NodeToDofGlue{VectorValue{3, Int32}}}, V::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Gridap.FESpaces.NodeToDofGlue{VectorValue{3, Int32}}}) @ Gridap.FESpaces C:\Users\zinli.julia\packages\Gridap\Qj36b\src\FESpaces\Assemblers.jl:285 [12] assemble_matrix(f::Function, U::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Gridap.FESpaces.NodeToDofGlue{VectorValue{3, Int32}}}, V::Gridap.FESpaces.UnconstrainedFESpace{Vector{Float64}, Gridap.FESpaces.NodeToDofGlue{VectorValue{3, Int32}}}) @ Gridap.FESpaces C:\Users\zinli.julia\packages\Gridap\Qj36b\src\FESpaces\Assemblers.jl:342 [13] top-level scope @ In[1]:26

[amartinhuertas]

Still the same error. So I guess the problem is with the outer, ⋅ or ∇ or ⊙

Not necessarily. As far I can see in the stack trace, the cell matrix which results from the evaluation of the bilinear form seems to be of the correct type, i.e., ComplexF64, while the data structure in which the matrix contributions are inserted during FE assembly does not (Gridap.Algebra.InserterCSC{Float64, Int64}). Let me investigate that further and I will get back to you.

[fverdugo]

@zlin-opt

You need to pass the key-word argument vector_type=Vector{ComplexF64} to the FESpace constructor

[zlin-opt]

Still the same error. So I guess the problem is with the outer, ⋅ or ∇ or ⊙

Not necessarily. As far I can see in the stack trace, the cell matrix which results from the evaluation of the bilinear form seems to be of the correct type, i.e., ComplexF64, while the data structure in which the matrix contributions are inserted during FE assembly does not (Gridap.Algebra.InserterCSC{Float64, Int64}). Let me investigate that further and I will get back to you.

thank

vector_type=Vector{ComplexF64}

great that works! I thought this line took care of making a complex FE space reffe = ReferenceFE(lagrangian,VectorValue{3,ComplexF64},order)

so I need to specify the type in the constructor itself

works now!

[amartinhuertas]

You need to pass the key-word argument vector_type=Vector{ComplexF64} to the FESpace constructor

Ok, perf. That was the cause of the type mismatch. In any case, AFAIK, mixing Floats and Complex numbers in the bilinear form can also cause trouble

[stevengj]

AFAIK, mixing Floats and Complex numbers in the bilinear form can also cause trouble

That hasn't been our experience IIRC — Julia's promotion rules should work fine for this.

But our PML tutorial did say:

Since our problem involves complex numbers (because of PML), we need to assign the vector_type to be Vector{ComplexF64}.

I thought this line took care of making a complex FE space reffe = ReferenceFE(lagrangian,VectorValue{3,ComplexF64},order)

I don't think you need or want the basis functions to be complex-valued. (Only the coefficients multiplying the basis functions are complex.) Doing this may make things more expensive than necessary.

[zlin-opt]

AFAIK, mixing Floats and Complex numbers in the bilinear form can also cause trouble

That hasn't been our experience IIRC — Julia's promotion rules should work fine for this.

But our PML tutorial did say:

Since our problem involves complex numbers (because of PML), we need to assign the vector_type to be Vector{ComplexF64}.

I thought this line took care of making a complex FE space reffe = ReferenceFE(lagrangian,VectorValue{3,ComplexF64},order)

I don't think you need or want the basis functions to be complex-valued. (Only the coefficients multiplying the basis functions are complex.) Doing this may make things more expensive than necessary.

I see. That makes sense.