Errors for FSI problem

[oriolcg]

Hi @santiagobadia, @ericneiva, @fverdugo

In this issue I'll report some of the errors that I'm facing when building a FSI driver using AgFEM.

In the following file you'll find the code that I'm executing.

https://gist.github.com/oriolcg/66f6fb257f77b674021671056a055039

It first builds the geometry of the elastic flag benchmark (simplified) and then I just set up an Stokes problem on the fluid domain. I'm linking with Gridap#master and GridapEmbedded#gridap_v0.16. In this preliminar setup I already have a couple of problems:

1. In line 123 I have the following call: Vᵥ_ser = FESpace(model_fluid,FEᵥ_ser)#,conformity=:L2). If I add the conformity, I get the following error

ERROR: LoadError: MethodError: no method matching Gridap.FESpaces.FESpace(::Gridap.Geometry.DiscreteModelPortion{2,2}, ::Gridap.FESpaces.CellFE{FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}}; constraint=nothing, conformity=:L2)
Closest candidates are:
  Gridap.FESpaces.FESpace(::Gridap.Geometry.DiscreteModel, ::Gridap.FESpaces.CellFE; labels, dirichlet_tags, dirichlet_masks, constraint, vector_type) at /Users/oriol/.julia/packages/Gridap/ONkqj/src/FESpaces/FESpaceFactories.jl:2 got unsupported keyword argument "conformity"
  Gridap.FESpaces.FESpace(::Gridap.Geometry.RestrictedDiscreteModel, ::Gridap.FESpaces.CellFE; constraint, kwargs...) at /Users/oriol/.julia/packages/Gridap/ONkqj/src/FESpaces/FESpaceFactories.jl:72
  Gridap.FESpaces.FESpace(::Gridap.Geometry.DiscreteModel, ::AbstractArray{var"#s351",N} where N where var"#s351"<:Gridap.ReferenceFEs.ReferenceFE; conformity, labels, dirichlet_tags, dirichlet_masks, constraint, vector_type) at /Users/oriol/.julia/packages/Gridap/ONkqj/src/FESpaces/FESpaceFactories.jl:83
  ...
Stacktrace:
 [1] kwerr(::NamedTuple{(:constraint, :conformity),Tuple{Nothing,Symbol}}, ::Type{T} where T, ::Gridap.Geometry.DiscreteModelPortion{2,2}, ::Gridap.FESpaces.CellFE{FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}}) at ./error.jl:157
 [2] Gridap.FESpaces.FESpace(::Gridap.Geometry.RestrictedDiscreteModel{2,2}, ::Gridap.FESpaces.CellFE{FillArrays.Fill{Int64,1,Tuple{Base.OneTo{Int64}}}}; constraint::Nothing, kwargs::Base.Iterators.Pairs{Symbol,Symbol,Tuple{Symbol},NamedTuple{(:conformity,),Tuple{Symbol}}}) at /Users/oriol/.julia/packages/Gridap/ONkqj/src/FESpaces/FESpaceFactories.jl:77
 [3] top-level scope at /Users/oriol/Progs/Gridap/GridapFSI.jl/tmp.jl:123
in expression starting at /Users/oriol/Progs/Gridap/GridapFSI.jl/tmp.jl:123

2. When adding the boundary terms in line 138: ∫( (γ/h)*v⋅u - v⋅(n_ΓDf⋅∇(u)) - (n_ΓDf⋅∇(v))⋅u + (p*n_ΓDf)⋅v + (q*n_ΓDf)⋅u )dΓDf, I get the following error:

ERROR: LoadError: MethodError: no method matching copy(::Gridap.Fields.ArrayBlock{Array{Float64,2},2})
Closest candidates are:
  copy(::TimerOutputs.TimeData) at /Users/oriol/.julia/packages/TimerOutputs/4QAIk/src/TimerOutput.jl:10
  copy(::Pkg.Types.VersionSpec) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.5/Pkg/src/versions.jl:230
  copy(::Revise.RelocatableExpr) at /Users/oriol/.julia/packages/Revise/qxX5H/src/relocatable_exprs.jl:25
  ...
Stacktrace:
 [1] _compress!(::Tuple{Gridap.Fields.ArrayBlock{Array{Float64,2},2},Gridap.Fields.ArrayBlock{Array{Float64,1},1}}, ::Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.CellData.ConstrainRowsMap,1,Tuple{Base.OneTo{Int64}}},Tuple{Gridap.Fields.ArrayBlock{Array{Float64,2},2},Gridap.Fields.ArrayBlock{Array{Float64,1},1}},1,Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.CellData.ConstrainColsMap,1,Tuple{Base.OneTo{Int64}}},Tuple{Gridap.Fields.ArrayBlock{Array{Float64,2},2},Gridap.Fields.ArrayBlock{Array{Float64,1},1}},1, ...

Any help will be very much appreciated

[ericneiva]

2. When adding the boundary terms in line 138: ∫( (γ/h)v⋅u - v⋅(n_ΓDf⋅∇(u)) - (n_ΓDf⋅∇(v))⋅u + (pn_ΓDf)⋅v + (q*n_ΓDf)⋅u )dΓDf, I get the following error:

For this one, we need to implement copy or _compress! for ArrayBlocks, isn't it?

[oriolcg]

Indeed, @ericneiva, there were some functions that had to be implemented for ArrayBlocks. But now I get the following error

ERROR: LoadError: MethodError: no method matching +(::Nothing, ::Nothing)
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:538
  +(::ChainRulesCore.DoesNotExist, ::Any) at /Users/oriol/.julia/packages/ChainRulesCore/D0go7/src/differential_arithmetic.jl:23
  +(::ChainRulesCore.Zero, ::Any) at /Users/oriol/.julia/packages/ChainRulesCore/D0go7/src/differential_arithmetic.jl:63
  ...
Stacktrace:
 [1] _compress!(::Tuple{Gridap.Fields.ArrayBlock{Array{Float64,2},2},Gridap.Fields.ArrayBlock{Array{Float64,1}

I'm not sure if this error comes from a bad implementation of the extended functions in the PR gridap/Gridap.jl#620 or for another reason in GridapEmbedded.

[ericneiva]

Hi, @oriolcg, I am taking a look at this now, because I have reproduced the error, while trying to let GridapEmbedded tests pass with the current Gridap#master...

[ericneiva]

I'm not sure if this error comes from a bad implementation of the extended functions in the PR gridap/Gridap.jl#620 or for another reason in GridapEmbedded.

I think you were trying to sum two array blocks in a way that it is not implemented. You can take a look at how I sum the array blocks, by inspecting the new _compress! implementation in gridap/Gridap.jl#623.

I have implemented setsize! for ArrayBlocks and _compress! for Tuple{ArrayBlock,ArrayBlock}. With this, I am able to solve again bimaterial poisson problems in GridapEmbedded.

But it is still not enough to pass all tests... I get the same error you report above @oriolcg with test/StokesAgFEMTests.jl.

In line 123 I have the following call: Vᵥ_ser = FESpace(model_fluid,FEᵥ_ser)#,conformity=:L2). If I add the conformity, I get the following error

[oriolcg]

Thanks @ericneiva, with the new implementation in gridap/Gridap.jl#623 I don't have the second error anymore.

Now I recover the error that I had some time ago:

ERROR: LoadError: AssertionError: A check failed
Stacktrace:
 [1] macro expansion at /Users/oriol/.julia/packages/Gridap/aDLdu/src/Helpers/Macros.jl:60 [inlined]
 [2] pos_and_neg_indices(::Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Reindex{Array{Int64,1}},1,Tuple{Base.OneTo{Int64}}},Int64,1,Tuple{Gridap.Arrays.LazyArray{FillArrays.Fill{Gridap.Arrays.Reindex{Gridap.Arrays.IdentityVector{Int64}},1,Tuple{Base.OneTo{Int64}}},Int64,1,Tuple{Array{Int64,1}}}}}) at /Users/oriol/.julia/packages/Gridap/aDLdu/src/Arrays/PosNegReindex.jl:112
 [3] lazy_map(::typeof(Gridap.Arrays.evaluate), ::Type{Array{Float64,2}}, ...

This error appears whenever there is at least one full inactive (solid) element within the flag.
This works:

This doesn't work:

Any idea on what could be the problem?

[amartinhuertas]

Any idea on what could be the problem?

I am not sure what is the particular cause of the problem in your case, but I have also faced this assertion error in other contexts.

The following post may help you: gridap/Gridap.jl#603 (comment)

[oriolcg]

Thanks @amartinhuertas. Then, it might be related to the fact that both sides of the interface are active when defining the AgFEMSpace. @santiagobadia, I think you mentioned that your PhD already solved a similar problem, did he/she face similar issues?

[fverdugo]

@oriolcg can you upload the mesh file so that I can try to reproduce the problem?

[oriolcg]

@fverdugo, the code is here: https://gist.github.com/oriolcg/66f6fb257f77b674021671056a055039
I'm using a CartesianDiscreteModel so no input mesh needed. The problem appears when the number of elements in the vertical direction is >12.

[fverdugo]

@oriolcg gridap/Gridap.jl#623 and gridap/Gridap.jl#624 fix the problem!

[oriolcg]

Thanks @fverdugo! I'll keep adding things and hope there are no more issues.

[fverdugo]

This should be fixed via the new release https://github.com/gridap/GridapEmbedded.jl/releases/tag/v0.7.0

If the error persists, reopen the issue.

[oriolcg]

Hi @fverdugo, @ericneiva, @santiagobadia,
I have a new issue with the FSI + AgFEM driver. When defining the FE space for the solid part with the same aggregates that I use for the fluid I get the following bounds error:

ERROR: LoadError: BoundsError: attempt to access 1170-element Array{Int64,1} at index [1171]
Stacktrace:
 [1] setindex! at ./array.jl:847 [inlined]
 [2] setindex! at ./multidimensional.jl:559 [inlined]
 [3] _prepare_DOF_to_DOFs(::Array{Int64,1}, ::Gridap.Arrays.Table{Int64,Array{Int64,1},Array{Int32,1}}, ::Gridap.Arrays.Table{Float64,Array{Float64,1},Array{Int32,1}}, ::Int64, ::Int64) at /Users/oriol/.julia/packages/Gridap/293eQ/src/FESpaces/FESpacesWithLinearConstraints.jl:114
 [4] Gridap.FESpaces.FESpaceWithLinearConstraints(::Array{Int64,1}, ::Gridap.Arrays.Table{Int64,Array{Int64,1},Array{Int32,1}}, ::Gridap.Arrays.Table{Float64,Array{Float64,1},Array{Int32,1}}, ::Gridap.FESpaces.ExtendedFESpace{Gridap.FESpaces.UnconstrainedFESpace{Array{Float64,1},Nothing}}) at /Users/oriol/.julia/packages/Gridap/293eQ/src/FESpaces/FESpacesWithLinearConstraints.jl:81
 [5] AgFEMSpace(::Gridap.FESpaces.ExtendedFESpace{Gridap.FESpaces.UnconstrainedFESpace{Array{Float64,1},Nothing}}, ::Array{Int32,1}, ::Gridap.FESpaces.SingleFieldFEBasis{Gridap.FESpaces.TestBasis,Gridap.CellData.PhysicalDomain}, ::Gridap.CellData.CellDof{Gridap.CellData.PhysicalDomain}) at /Users/oriol/.julia/packages/GridapEmbedded/mBCb5/src/AgFEM/AgFEMSpaces.jl:62
 [6] AgFEMSpace(::Gridap.FESpaces.ExtendedFESpace{Gridap.FESpaces.UnconstrainedFESpace{Array{Float64,1},Nothing}}, ::Array{Int32,1}, ::Gridap.FESpaces.ExtendedFESpace{Gridap.FESpaces.UnconstrainedFESpace{Array{Float64,1},Nothing}}) at /Users/oriol/.julia/packages/GridapEmbedded/mBCb5/src/AgFEM/AgFEMSpaces.jl:3
 [7] top-level scope at /Users/oriol/Progs/Gridap/GridapFSI.jl/tmp.jl:142
in expression starting at /Users/oriol/Progs/Gridap/GridapFSI.jl/tmp.jl:142

The error appears in line 142 of the following code: https://gist.github.com/oriolcg/636441a6d9e2420165033f50dc0f15f3#file-agfem_fsi-jl-L142.

Do I need to define other aggregates? or a different aggregation strategy?

[fverdugo]

You need different aggregates in each case. Fluid and solid fe spaces are defined over different domains.

[ericneiva]

This incomplete snippet of a bimaterial Poisson problem might help you:

geom⁺ = geom; geom⁻ = ! geom⁺
cutgeo = cut(bgmodel,union(geom⁺,geom⁻))
. . .
aggs⁺ = aggregate(strategy,cutgeo,geom⁺)
aggs⁻ = aggregate(strategy,cutgeo,geom⁻)
. . . 
V⁺ = AgFEMSpace(Vstd⁺,aggs⁺)
V⁻ = AgFEMSpace(Vstd⁻,aggs⁻)