Skip to content

Commit

Permalink
More fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
JordiManyer committed Aug 14, 2023
1 parent 9b1620c commit 265a9d0
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 34 deletions.
2 changes: 1 addition & 1 deletion src/FESpaces.jl
Original file line number Diff line number Diff line change
Expand Up @@ -489,7 +489,7 @@ function _find_vector_type(spaces,gids)
# TODO Now the user can select the local vector type but not the global one
# new kw-arg global_vector_type ?
# we use PVector for the moment
local_vector_type = get_vector_type(PartitonedArrays.getany(spaces))
local_vector_type = get_vector_type(PartitionedArrays.getany(spaces))

if local_vector_type <: BlockVector
T = eltype(local_vector_type)
Expand Down
68 changes: 37 additions & 31 deletions src/MultiField.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,11 +37,7 @@ struct DistributedMultiFieldFESpace{MS,A,B,C,D} <: DistributedFESpace
vector_type::Type{D}
function DistributedMultiFieldFESpace(
field_fe_space::AbstractVector{<:DistributedSingleFieldFESpace},
<<<<<<< HEAD
part_fe_space::AbstractPData{<:MultiFieldFESpace{MS}},
=======
part_fe_space::AbstractArray{<:MultiFieldFESpace},
>>>>>>> 38b02ee2811f5a28d66359dd085b826041f9ee07
part_fe_space::AbstractArray{<:MultiFieldFESpace{MS}},
gids::PRange,
vector_type::Type{D}) where {D,MS}
A = typeof(field_fe_space)
Expand Down Expand Up @@ -425,52 +421,59 @@ end

function FESpaces.SparseMatrixAssembler(
trial::DistributedMultiFieldFESpace{<:BlockMultiFieldStyle},
test::DistributedMultiFieldFESpace{<:BlockMultiFieldStyle},
test ::DistributedMultiFieldFESpace{<:BlockMultiFieldStyle},
par_strategy=SubAssembledRows())
Tv = get_vector_type(get_part(local_views(first(trial))))
Tv = get_vector_type(PartitionedArrays.getany(local_views(first(trial))))
T = eltype(Tv)
Tm = SparseMatrixCSC{T,Int}
SparseMatrixAssembler(Tm,Tv,trial,test,par_strategy)
end

function local_views(a::MultiField.BlockSparseMatrixAssembler{NB,NV,SB,P}) where {NB,NV,SB,P}
assems = a.block_assemblers
parts = get_part_ids(local_views(first(assems)))
map(parts) do p
idx = CartesianIndices(axes(assems))
block_assems = map(idx) do I
get_part(local_views(assems[I]),p)
# Array of PArrays -> PArray of Arrays
function to_parray_of_arrays(a::AbstractArray{<:MPIArray})
indices = linear_indices(first(a))
map(indices) do i
map(a) do aj
PartitionedArrays.getany(aj)
end
return MultiField.BlockSparseMatrixAssembler{NB,NV,SB,P}(block_assems)
end
end

function local_views(a::MatrixBlock,rows,cols)
parts = get_part_ids(local_views(first(a.array)))
map(parts) do p
idx = CartesianIndices(axes(a))
array = map(idx) do I
get_part(local_views(a[I],rows[I[1]],cols[I[2]]),p)
function to_parray_of_arrays(a::AbstractArray{<:DebugArray})
indices = linear_indices(first(a))
map(indices) do i
map(a) do aj
aj.items[i]
end
ArrayBlock(array,a.touched)
end
end

function local_views(a::MultiField.BlockSparseMatrixAssembler{NB,NV,SB,P}) where {NB,NV,SB,P}
assems = a.block_assemblers
array = to_parray_of_arrays(map(local_views,assems))
return map(MultiField.BlockSparseMatrixAssembler{NB,NV,SB,P},array)
end

function local_views(a::MatrixBlock,rows,cols)
idx = CartesianIndices(axes(a))
array = map(idx) do I
local_views(a[I],rows[I[1]],cols[I[2]])
end
return map(b -> ArrayBlock(b,a.touched), to_parray_of_arrays(array))
end

function local_views(a::VectorBlock,rows)
parts = get_part_ids(local_views(first(a.array)))
map(parts) do p
idx = CartesianIndices(axes(a))
array = map(idx) do I
get_part(local_views(a[I],rows[I]),p)
end
ArrayBlock(array,a.touched)
idx = CartesianIndices(axes(a))
array = map(idx) do I
local_views(a[I],rows[I])
end
return map(b -> ArrayBlock(b,a.touched), to_parray_of_arrays(array))
end

function local_views(a::ArrayBlockView,axes...)
array = local_views(a.array,axes...)
map(array) do array
MultiField.ArrayBlockView(array,a.block_map)
ArrayBlockView(array,a.block_map)
end
end

Expand Down Expand Up @@ -514,10 +517,12 @@ end
#! but requires a little bit of refactoring of the assembly code. Postponed until GridapDistributed v0.3.

function Algebra.create_from_nz(a::ArrayBlock{<:DistributedAllocationCOO{<:FullyAssembledRows}})
array = map(_fa_create_from_nz_temporary_fix,a.array)
#array = map(_fa_create_from_nz_temporary_fix,a.array)
array = map(Algebra.create_from_nz,a.array)
return mortar(array)
end

"""
function _fa_create_from_nz_temporary_fix(a::DistributedAllocationCOO{<:FullyAssembledRows})
parts = get_part_ids(local_views(a))
Expand Down Expand Up @@ -565,3 +570,4 @@ function _fa_create_from_nz_temporary_fix(a::DistributedAllocationCOO{<:FullyAss
exchanger = empty_exchanger(parts)
return PSparseMatrix(values,rows,cols,exchanger)
end
"""
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
module BlockSparseMatrixAssemblersTests

using Test, LinearAlgebra, BlockArrays, SparseArrays

using Gridap
Expand All @@ -6,11 +8,14 @@ using Gridap.FESpaces, Gridap.ReferenceFEs, Gridap.MultiField
using GridapDistributed
using PartitionedArrays

parts = get_part_ids(SequentialBackend(),(2,2))
nparts = (2,2)
parts = with_debug() do distribute
distribute(LinearIndices((prod(nparts),)))
end

sol(x) = sum(x)

model = CartesianDiscreteModel(parts,(0.0,1.0,0.0,1.0),(12,12))
model = CartesianDiscreteModel(parts,nparts,(0.0,1.0,0.0,1.0),(12,12))
Ω = Triangulation(model)

reffe = LagrangianRefFE(Float64,QUAD,1)
Expand Down Expand Up @@ -106,6 +111,11 @@ block_trials = map(range -> get_block_fespace(X.field_fe_space,range),block_rang

assem_blocks = SparseMatrixAssembler(Xb,Yb,FullyAssembledRows())

local_views(assem_blocks)

ab = assem_blocks.block_assemblers
map(local_views,ab)

A1_blocks = assemble_matrix(assem_blocks,bmatdata);
b1_blocks = assemble_vector(assem_blocks,bvecdata);

Expand All @@ -131,3 +141,4 @@ A11 = A1_blocks.blocks[1,1]
A12 = A1_blocks.blocks[1,2]
A22 = A1_blocks.blocks[2,2]

end

0 comments on commit 265a9d0

Please sign in to comment.