Issue with DG-SIP and FullyAssembledRows

[amartinhuertas]

When activating FullyAssembledRows assembly strategy in the DG-SIP method, e.g., in PoissonTests.jl, the code breaks during FE assembly (i.e., AffineFEOperator constructor)

The code below reproduces the issue:

module PoissonTests
using SparseMatricesCSR
using Gridap
using Gridap.Algebra
using Gridap.FESpaces
using GridapDistributed
using PartitionedArrays
using Test

function main(distribute, parts)
  ranks = distribute(LinearIndices((prod(parts),)))
  output = mkpath(joinpath(@__DIR__,"output"))

  domain = (0,4,0,4)
  cells = (4,4)
  model = CartesianDiscreteModel(ranks, parts,domain,cells)

  labels = get_face_labeling(model)
  add_tag_from_tags!(labels,"dirichlet",[1,2,3,5,7])
  add_tag_from_tags!(labels,"neumann",[4,6,8])

  Ω = Triangulation(with_ghost,model)
  Γn = Boundary(with_ghost,model,tags="neumann")
  n_Γn = get_normal_vector(Γn)

  k = 2
  u((x,y)) = (x+y)^k
  f(x) = -Δ(u,x)
  g = n_Γn⋅∇(u)

  reffe = ReferenceFE(lagrangian,Float64,k)
  V = TestFESpace(model,reffe,dirichlet_tags="dirichlet")
  U = TrialFESpace(u,V)

  dΩ = Measure(Ω,2*k)
  dΓn = Measure(Γn,2*k)

  a(u,v) = ∫( ∇(v)⋅∇(u) )dΩ
  l(v) = ∫( v*f )dΩ + ∫( v*g )dΓn
  assem=SparseMatrixAssembler(SparseMatrixCSR{0,Float64,Int},Vector{Float64},U,V,FullyAssembledRows())
  op = AffineFEOperator(a,l,U,V,assem)

  uh = solve(op)
  eh = u - uh
  @test sqrt(sum( ∫(abs2(eh))dΩ )) < 1.0e-9

  # Now with DG
  h = 2
  γ = 10

  V_dg = FESpace(model,reffe,conformity=:L2)
  U_dg = TrialFESpace(V_dg)

  Λ = Skeleton(with_ghost,model)
  Γd = Boundary(with_ghost,model,tags="dirichlet")

  dΛ = Measure(Λ,2*k)
  dΓd = Measure(Γd,2*k)

  n_Γd = get_normal_vector(Γd)
  n_Λ = get_normal_vector(Λ)

  a_dg(u,v) =
    ∫( ∇(v)⋅∇(u) )*dΩ +
    ∫( (γ/h)*v*u  - v*(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))*u )*dΓd +
    ∫( (γ/h)*jump(v*n_Λ)⋅jump(u*n_Λ) -
       jump(v*n_Λ)⋅mean(∇(u)) -
       mean(∇(v))⋅jump(u*n_Λ) )*dΛ

  l_dg(v) =
    ∫( v*f )*dΩ +
    ∫( v*g )dΓn +
    ∫( (γ/h)*v*u - (n_Γd⋅∇(v))*u )*dΓd


  assem=SparseMatrixAssembler(U_dg,V_dg, FullyAssembledRows())
  op = AffineFEOperator(a_dg,l_dg,U_dg,V_dg,assem)
  uh = solve(op)
  eh = u - uh
  @test sqrt(sum( ∫(abs2(eh))dΩ )) < 1.0e-9

end
end # module

when executed, e.g., as:

module PoissonTestsSeq
using PartitionedArrays
include("../PoissonTests.jl")
with_debug() do distribute
    PoissonTests.main(distribute,(2,2))
end 
end # module

[amartinhuertas]

@JordiManyer ... where you aware of this?

[JordiManyer]

This might be because off a recent change. In any case, I don't thing this worked anyways for the FullyAssembledRows assembly, since is has jumps (and therefore will always get non-local contributions).

[amartinhuertas]

Hi @JordiManyer ... thanks for your answer. Let me answer to you below.

In any case, I don't thing this worked anyways for the FullyAssembledRows assembly, since is has jumps (and therefore will always get non-local contributions).

I agree that FullyAssembledRows does not work with DG-SIP and, e.g., a H-div conforming discretization (as we discussed a while ago). However, if I am not wrong, it should work with piece-wise discontinuous polynomials, as in the tutorial.

Anyways, I believe I have been able to narrow down the problem to a MWE. The following Gridap snippet:

using Gridap
using FillArrays
model = CartesianDiscreteModel((0,1,0,1),(3,3))
trian = BoundaryTriangulation(model,Int[])
reffe = ReferenceFE(lagrangian,Float64,2)
V  = FESpace(model,reffe)
uh = zero(V)
cellvals = get_cell_dof_values(uh,trian)
println(length(Gridap.Arrays.testitem(cellvals)))

weirdly generates an array of size 9. This causes the error in the GridapDistributed.jl script above, in particular, for the particular processor/part that is not in contact with the Dirichlet boundary. For completeness, the stack trace for the error in the script is attached below. Note the number 9 in the ERROR message in the vector x. This vector exactly comes from the last call in the mwe above.

Does it ring a bell? Have you ever seen an error like this? Thanks!

ERROR: DimensionMismatch: second dimension of A, 1, does not match length of x, 9
Stacktrace:
  [1] gemv!(y::Vector{Float64}, tA::Char, A::Matrix{Float64}, x::Vector{Float64}, α::Bool, β::Bool)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:404
  [2] generic_matvecmul!
    @ ~/.julia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:71 [inlined]
  [3] mul!
    @ ~/.julia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:66 [inlined]
  [4] mul!
    @ ~/.julia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:237 [inlined]
  [5] *
    @ ~/.julia/juliaup/julia-1.10.10+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/matmul.jl:53 [inlined]
  [6] return_cache(k::Gridap.Arrays.MulAddMap{Int64}, a::Matrix{Float64}, b::Vector{Float64}, c::Vector{Float64})
    @ Gridap.Arrays ~/.julia/packages/Gridap/6UkFZ/src/Arrays/AlgebraMaps.jl:65
  [7] evaluate
    @ ~/.julia/packages/Gridap/6UkFZ/src/Arrays/Maps.jl:86 [inlined]
  [8] return_value
    @ ~/.julia/packages/Gridap/6UkFZ/src/Arrays/Maps.jl:64 [inlined]
  [9] return_value
    @ ~/.julia/packages/Gridap/6UkFZ/src/CellData/AttachDirichlet.jl:64 [inlined]
 [10] return_type
    @ ~/.julia/packages/Gridap/6UkFZ/src/Arrays/Maps.jl:62 [inlined]
 [11] lazy_map(::Gridap.CellData.AttachDirichletMap, ::Gridap.Arrays.LazyArray{…}, ::Gridap.Arrays.LazyArray{…}, ::Fill{…})
    @ Gridap.Arrays ~/.julia/packages/Gridap/6UkFZ/src/Arrays/LazyArrays.jl:57
 [12] attach_dirichlet(cellmatvec::Gridap.Arrays.LazyArray{…}, cellvals::Gridap.Arrays.LazyArray{…}, cellmask::Fill{…})
    @ Gridap.CellData ~/.julia/packages/Gridap/6UkFZ/src/CellData/AttachDirichlet.jl:4
 [13] _pair_contribution_when_possible(biform::Gridap.CellData.DomainContribution, liform::Gridap.CellData.DomainContribution, uhd::Gridap.FESpaces.SingleFieldFEFunction{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/6UkFZ/src/FESpaces/Assemblers.jl:515
 [14] collect_cell_matrix_and_vector
    @ ~/.julia/packages/Gridap/6UkFZ/src/FESpaces/Assemblers.jl:482 [inlined]
 [15] #4
    @ ./generator.jl:36 [inlined]
 [16] iterate
    @ ./generator.jl:47 [inlined]
 [17] collect_to!(dest::Vector{…}, itr::Base.Generator{…}, offs::Int64, st::NTuple{…})
    @ Base ./array.jl:892
 [18] collect_to_with_first!
    @ ./array.jl:870 [inlined]
 [19] collect(itr::Base.Generator{Base.Iterators.Zip{…}, Base.var"#4#5"{…}})
    @ Base ./array.jl:844
 [20] map
    @ ./abstractarray.jl:3417 [inlined]
 [21] map
    @ ~/.julia/packages/PartitionedArrays/MVmxR/src/debug_array.jl:105 [inlined]
 [22] collect_cell_matrix_and_vector(trial::GridapDistributed.DistributedSingleFieldFESpace{…}, test::GridapDistributed.DistributedSingleFieldFESpace{…}, biform::GridapDistributed.DistributedDomainContribution{…}, liform::GridapDistributed.DistributedDomainContribution{…}, uhd::GridapDistributed.DistributedCellField{…})
    @ GridapDistributed ~/git-repos/GridapDistributed.jl/src/FESpaces.jl:617
 [23] AffineFEOperator(weakform::Gridap.FESpaces.var"#84#85"{…}, trial::GridapDistributed.DistributedSingleFieldFESpace{…}, test::GridapDistributed.DistributedSingleFieldFESpace{…}, assem::GridapDistributed.DistributedSparseMatrixAssembler{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/6UkFZ/src/FESpaces/AffineFEOperators.jl:37
 [24] AffineFEOperator(::Main.PoissonTestsSeq.PoissonTests.var"#a_dg#5"{…}, ::Main.PoissonTestsSeq.PoissonTests.var"#l_dg#6"{…}, ::GridapDistributed.DistributedSingleFieldFESpace{…}, ::Vararg{…})
    @ Gridap.FESpaces ~/.julia/packages/Gridap/6UkFZ/src/FESpaces/AffineFEOperators.jl:51
 [25] main(distribute::Type{…}, parts::Tuple{…}, parallel_assembly_strategy::GridapDistributed.FullyAssembledRows)
    @ Main.PoissonTestsSeq.PoissonTests ~/git-repos/GridapDistributed.jl/test/PoissonTests.jl:92
 [26] #1
    @ ~/git-repos/GridapDistributed.jl/test/sequential/PoissonTests.jl:6 [inlined]
 [27] with_debug(f::Main.PoissonTestsSeq.var"#1#2")
    @ PartitionedArrays ~/.julia/packages/PartitionedArrays/MVmxR/src/debug_array.jl:8
 [28] top-level scope
    @ ~/git-repos/GridapDistributed.jl/test/sequential/PoissonTests.jl:4
Some type information was truncated. Use `show(err)` to see complete types.

[JordiManyer]

The change I did is the following: #180 the other section.

The idea is that before, BoundaryTriangulation was returning the local processor boundary, not the global boundary. I changed it to returning the global boundary portion within the processor. Indeed this creates empty triangulations in fully interior processors.

Despite this, I think it was the correct thing to do. It is not intuitive (and I dont think its what anyone expects) to obtain the local processor boundary like we had before. It also simplifies and optimizes the whole process of adding the ghosts back. Both Zach and Santi agreed on this.

The issue, I guess, is that now there is no constructor for getting the local processor boundary which is what you now need for SIP. You want to capture half of the jumps on the processor boundary, so you need to be integrating there. But it should not be called a BoundaryTriangulation.

Basically what you want is this:

trians = map(BoundaryTriangulation,local_views(model)) 
DistributedTriangulation(trians, model)

In any case, I would just not use FullyAssembledRows.

[amartinhuertas]

The change I did is the following: #180 the other section.

Are you 100% sure there is cause-and-effect among this change and the error?

Despite this, I think it was the correct thing to do. It is not intuitive (and I dont think its what anyone expects) to obtain the local processor boundary like we had before. It also simplifies and optimizes the whole process of adding the ghosts back. Both Zach and Santi agreed on this.

Ok. But I am not sure if this applies to this particular case. We are not using Boundary(model), but Boundary(model, tags="dirichlet"), and I believe that the tags have always been respecting that of the global model.

The issue, I guess, is that now there is no constructor for getting the local processor boundary which is what you now need for SIP. You want to capture half of the jumps on the processor boundary, so you need to be integrating there. But it should not be called a BoundaryTriangulation.

I don't think I need that (the local processor boundary).

I need, for the skeleton triangulation, all skeleton facets that either:

• (1) are shared among owned cells.
• (2) are shared among owned and ghost cells, or the other way round.
• (3) are shared among ghost cells.

Facets that are only on ghost cells should not be on the local skeleton triangulation.

I believe the code is already fullfilling this.

Why I do not current get, is why we are experiencing the error with FullyAssembledRows, and not with SubAssembledRows.

[amartinhuertas]

And finally, I think it is an issue with Gridap related to the way we handle zero-size arrays, is not a GridapDistributed issue.

[JordiManyer]

We are not using Boundary(model), but Boundary(model, tags="dirichlet"), and I believe that the tags have always been respecting that of the global model.

Ok, then yeah there is no change.

[amartinhuertas]

A follow up on this one. It turns out that the difference among SubAssembledRows and FullyAssembledRows is in the t array passed to the attach_dirichlet call in the _pair_contribution_when_possible(biform,liform,uhd) function:

matvec.dict[trian] = attach_dirichlet(t,cellvals,cellmask)

With SubAssembledRows, tesitem(t) returns arrays of compatible shape and dimensions to cellvals. However, with FullyAssembledRows, it does not. Thus, as apposed to what I said in a post before, it seems that the size of the cellvals array being 9 seems to be appropriate.

On the other hand, I have been able to simplify/narrow down the bilinear and linear forms in the code corresponding to the first post to the following minimal version that still reproduces the error:

a_dg(u,v) =
    ∫( ∇(v)⋅∇(u) )*dΩ +
    ∫( v*u )*dΓd

l_dg(v) =
    ∫( v*f )*dΩ +
    ∫( v*u )*dΓd

[JordiManyer]

@amartinhuertas What version of Gridap are you using? For me, this

using Gridap
using FillArrays
model = CartesianDiscreteModel((0,1,0,1),(3,3))
trian = BoundaryTriangulation(model,Int[])
reffe = ReferenceFE(lagrangian,Float64,2)
V  = FESpace(model,reffe)
uh = zero(V)
cellvals = get_cell_dof_values(uh,trian)
println(length(Gridap.Arrays.testitem(cellvals)))

returns an array of size 0, as expected.

So it might be that you have broken something in your own branch, or that you are using an old version of Gridap?
FYI: I am testing from the master branch.

[amartinhuertas]

Hi @JordiManyer,

that's interesting.

I am using Gridap 0.19.6 (without touching anything I am aware of). I indeed removed it from the package manager, and add it again, and same behaviour, i.e., an array of size 9.

I also tried just now with the latest version as of the master branch (9cf59d78cd79b7ededc43a73a8393131118cc045 from 7/11), and I also get an array of size 9.

So I am not able to reproduce your behaviour.

[amartinhuertas]

btw ... I am using julia version 1.10.10, not sure if this is the source of any difference.

As of the rest of direct dependencies of Gridap:

Status `~/git-repos/Gridap.jl/Project.toml`
  [1520ce14] AbstractTrees v0.4.5
  [fbb218c0] BSON v0.3.9
⌃ [8e7c35d0] BlockArrays v1.7.2
  [861a8166] Combinatorics v1.0.3
⌃ [864edb3b] DataStructures v0.19.1
  [ffbed154] DocStringExtensions v0.9.5
⌃ [442a2c76] FastGaussQuadrature v1.0.2
⌃ [5789e2e9] FileIO v1.17.0
⌃ [1a297f60] FillArrays v1.13.0
⌃ [f6369f11] ForwardDiff v1.1.0
⌅ [033835bb] JLD2 v0.5.15
⌅ [682c06a0] JSON v0.21.4
  [d3d80556] LineSearches v7.4.0
  [2774e3e8] NLsolve v4.5.1
  [b8a86587] NearestNeighbors v0.4.22
  [c74db56a] PolynomialBases v0.4.25
  [21216c6a] Preferences v1.5.0
  [1fd47b50] QuadGK v2.11.2
  [a0a7dd2c] SparseMatricesCSR v0.6.9
  [90137ffa] StaticArrays v1.9.15
  [64499a7a] WriteVTK v1.21.2
  [37e2e46d] LinearAlgebra
  [9a3f8284] Random
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0
  [8dfed614] Test

[amartinhuertas]

I have tried to update all Gridap dependencies to see if this makes any different, but it does NOT

Project Gridap v0.19.6
Status `~/git-repos/Gridap.jl/Project.toml`
  [1520ce14] AbstractTrees v0.4.5
  [fbb218c0] BSON v0.3.9
  [8e7c35d0] BlockArrays v1.9.1
  [861a8166] Combinatorics v1.0.3
  [864edb3b] DataStructures v0.19.3
  [ffbed154] DocStringExtensions v0.9.5
  [442a2c76] FastGaussQuadrature v1.1.0
  [5789e2e9] FileIO v1.17.1
  [1a297f60] FillArrays v1.15.0
  [f6369f11] ForwardDiff v1.2.2
⌅ [033835bb] JLD2 v0.5.15
⌅ [682c06a0] JSON v0.21.4
  [d3d80556] LineSearches v7.4.0
  [2774e3e8] NLsolve v4.5.1
  [b8a86587] NearestNeighbors v0.4.22
  [c74db56a] PolynomialBases v0.4.25
  [21216c6a] Preferences v1.5.0
  [1fd47b50] QuadGK v2.11.2
  [a0a7dd2c] SparseMatricesCSR v0.6.9
  [90137ffa] StaticArrays v1.9.15
  [64499a7a] WriteVTK v1.21.2
  [37e2e46d] LinearAlgebra
  [9a3f8284] Random
  [2f01184e] SparseArrays v1.10.0
  [10745b16] Statistics v1.10.0
  [8dfed614] Test
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

[JordiManyer]

Ok, sorry. My bad. I was not checking testitem, just the length of the array.

[amartinhuertas]

No worries. At least we have the same behaviour now. I am not sure if the length of that array is wrong. Indeed, it is returning the same thing for SubAssembledRows, that works flawlessly. As I said above, I believe is the array that comes from pair_contribution_if_possible in the case of FullyAssembledRows. I will keep investigating this.