From 14b54271646c78876397ce634bc72148096d0997 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 11 Dec 2023 12:38:12 +1100 Subject: [PATCH 01/25] Started porting code to new version of Gridap.ODEs --- src/GridapDistributed.jl | 4 +--- src/TransientDistributedCellField.jl | 10 +++++----- src/TransientFESpaces.jl | 20 +++++++++---------- ...TransientMultiFieldDistributedCellField.jl | 4 ++-- test/StokesOpenBoundaryTests.jl | 8 ++++++-- test/TransientDistributedCellFieldTests.jl | 4 ++-- ...ientMultiFieldDistributedCellFieldTests.jl | 6 +++--- 7 files changed, 29 insertions(+), 27 deletions(-) diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl index 27f4fc02..93e1077c 100644 --- a/src/GridapDistributed.jl +++ b/src/GridapDistributed.jl @@ -12,8 +12,7 @@ using Gridap.CellData using Gridap.Visualization using Gridap.FESpaces using Gridap.MultiField -using Gridap.ODEs.TransientFETools -using Gridap.ODEs.ODETools +using Gridap.ODEs using MPI using PartitionedArrays @@ -29,7 +28,6 @@ import Gridap.TensorValues: inner, outer, double_contraction, symmetric_part import LinearAlgebra: det, tr, cross, dot, ⋅, diag import Base: inv, abs, abs2, *, +, -, /, adjoint, transpose, real, imag, conj, getproperty, propertynames import Gridap.Fields: grad2curl -import Gridap.ODEs.ODETools: ∂t, ∂tt export FullyAssembledRows export SubAssembledRows diff --git a/src/TransientDistributedCellField.jl b/src/TransientDistributedCellField.jl index 8934f40f..3186d986 100644 --- a/src/TransientDistributedCellField.jl +++ b/src/TransientDistributedCellField.jl @@ -10,21 +10,21 @@ end local_views(f::TransientSingleFieldDistributedCellField) = local_views(f.cellfield) # Constructors -function TransientFETools.TransientCellField(single_field::DistributedSingleFieldFEFunction,derivatives::Tuple) +function ODEs.TransientCellField(single_field::DistributedSingleFieldFEFunction,derivatives::Tuple) TransientSingleFieldDistributedCellField(single_field,derivatives) end -function TransientFETools.TransientCellField(single_field::DistributedCellField,derivatives::Tuple) -TransientSingleFieldDistributedCellField(single_field,derivatives) +function ODEs.TransientCellField(single_field::DistributedCellField,derivatives::Tuple) + TransientSingleFieldDistributedCellField(single_field,derivatives) end # Time derivative -function ∂t(f::TransientDistributedCellField) +function ODEs.∂t(f::TransientDistributedCellField) cellfield, derivatives = first_and_tail(f.derivatives) TransientCellField(cellfield,derivatives) end -∂tt(f::TransientDistributedCellField) = ∂t(∂t(f)) +ODEs.∂tt(f::TransientDistributedCellField) = ∂t(∂t(f)) # Integration related function Fields.integrate(f::TransientDistributedCellField,b::DistributedMeasure) diff --git a/src/TransientFESpaces.jl b/src/TransientFESpaces.jl index 938d86c8..6b57a014 100644 --- a/src/TransientFESpaces.jl +++ b/src/TransientFESpaces.jl @@ -4,10 +4,10 @@ Fields.evaluate(U::DistributedSingleFieldFESpace,t::Nothing) = U (U::DistributedSingleFieldFESpace)(t) = U -∂t(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U) -∂t(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂t.(U.field_fe_space)) -∂tt(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U) -∂tt(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.field_fe_spaces)) +ODEs.∂t(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U) +ODEs.∂t(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂t.(U.field_fe_space)) +ODEs.∂tt(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U) +ODEs.∂tt(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.field_fe_spaces)) function TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace}) MultiFieldFESpace(spaces) @@ -15,24 +15,24 @@ end # Functions for transient FE Functions -function ODETools.allocate_jacobian( - op::TransientFETools.TransientFEOperatorFromWeakForm, +function ODEs.allocate_jacobian( + op::ODEs.TransientFEOpFromWeakForm, t0::Real, duh::Union{DistributedCellField,DistributedMultiFieldFEFunction}, cache) - _matdata_jacobians = TransientFETools.fill_initial_jacobians(op,t0,duh) + _matdata_jacobians = ODEs.fill_initial_jacobians(op,t0,duh) matdata = _vcat_distributed_matdata(_matdata_jacobians) allocate_matrix(op.assem_t,matdata) end -function ODETools.jacobians!( +function ODEs.jacobians!( A::AbstractMatrix, - op::TransientFETools.TransientFEOperatorFromWeakForm, + op::ODEs.TransientFEOpFromWeakForm, t::Real, xh::TransientDistributedCellField, γ::Tuple{Vararg{Real}}, cache) - _matdata_jacobians = TransientFETools.fill_jacobians(op,t,xh,γ) + _matdata_jacobians = ODEs.fill_jacobians(op,t,xh,γ) matdata = _vcat_distributed_matdata(_matdata_jacobians) assemble_matrix_add!(A,op.assem_t, matdata) A diff --git a/src/TransientMultiFieldDistributedCellField.jl b/src/TransientMultiFieldDistributedCellField.jl index 4bd280db..3f4fd47f 100644 --- a/src/TransientMultiFieldDistributedCellField.jl +++ b/src/TransientMultiFieldDistributedCellField.jl @@ -8,7 +8,7 @@ end local_views(f::TransientMultiFieldDistributedCellField) = local_views(f.cellfield) # Constructors -function TransientFETools.TransientCellField(multi_field::DistributedMultiFieldFEFunction,derivatives::Tuple) +function ODEs.TransientCellField(multi_field::DistributedMultiFieldFEFunction,derivatives::Tuple) transient_single_fields = _to_transient_single_distributed_fields(multi_field,derivatives) TransientMultiFieldDistributedCellField(multi_field,derivatives,transient_single_fields) end @@ -53,7 +53,7 @@ Base.iterate(f::TransientMultiFieldDistributedCellField) = iterate(f.transient_ Base.iterate(f::TransientMultiFieldDistributedCellField,state) = iterate(f.transient_single_fields,state) # Time derivative -function ∂t(f::TransientMultiFieldDistributedCellField) +function ODEs.∂t(f::TransientMultiFieldDistributedCellField) cellfield, derivatives = first_and_tail(f.derivatives) transient_single_field_derivatives = TransientDistributedCellField[] for transient_single_field in f.transient_single_fields diff --git a/test/StokesOpenBoundaryTests.jl b/test/StokesOpenBoundaryTests.jl index 8431e89b..93bf733a 100644 --- a/test/StokesOpenBoundaryTests.jl +++ b/test/StokesOpenBoundaryTests.jl @@ -81,7 +81,7 @@ function main(distribute,parts) ph0 = interpolate_everywhere(p(0.0),P0) xh0 = interpolate_everywhere([uh0,ph0],X0) - op = TransientFEOperator(res,jac,jac_t,X,Y) + op = TransientFEOperator(res,(jac,jac_t),X,Y) t0 = 0.0 tF = 1.0 @@ -90,7 +90,7 @@ function main(distribute,parts) ls = LUSolver() ode_solver = ThetaMethod(ls,dt,θ) - sol_t = solve(ode_solver,op,xh0,t0,tF) + sol_t = solve(ode_solver,op,t0,tF,xh0) l2(w) = w⋅w @@ -112,4 +112,8 @@ function main(distribute,parts) end end +with_debug() do distribute + main(distribute,(2,2)) +end + end #module diff --git a/test/TransientDistributedCellFieldTests.jl b/test/TransientDistributedCellFieldTests.jl index e356abe4..aed283d1 100644 --- a/test/TransientDistributedCellFieldTests.jl +++ b/test/TransientDistributedCellFieldTests.jl @@ -2,8 +2,8 @@ module TransientDistributedCellFieldTests using Gridap using GridapDistributed -using Gridap.ODEs.ODETools: ∂t, ∂tt -using Gridap.ODEs.TransientFETools: TransientCellField +using Gridap.ODEs: ∂t, ∂tt +using Gridap.ODEs: TransientCellField using PartitionedArrays using Test diff --git a/test/TransientMultiFieldDistributedCellFieldTests.jl b/test/TransientMultiFieldDistributedCellFieldTests.jl index bdb6eca6..76b66a76 100644 --- a/test/TransientMultiFieldDistributedCellFieldTests.jl +++ b/test/TransientMultiFieldDistributedCellFieldTests.jl @@ -2,9 +2,9 @@ module TransientMultiFieldDistributedCellFieldTests using Gridap using GridapDistributed -using Gridap.ODEs.ODETools: ∂t, ∂tt -using Gridap.ODEs.TransientFETools: TransientCellField -using Gridap.ODEs.TransientFETools: TransientTrialFESpace, TransientMultiFieldFESpace +using Gridap.ODEs: ∂t, ∂tt +using Gridap.ODEs: TransientCellField +using Gridap.ODEs: TransientTrialFESpace, TransientMultiFieldFESpace using PartitionedArrays using Test From c55d1aa3b53758ada17bdcc8b2f20ff76919bf75 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 19 Dec 2023 14:25:39 +1100 Subject: [PATCH 02/25] DistributedCellField now inherits from CellField To accomodate for this changes, we now keep a pointer to the DistributedTriangulation throughout the FESpace, CellField and Cellpoint. --- src/CellData.jl | 205 ++++++++++++++++----------- src/FESpaces.jl | 48 ++++--- src/GridapDistributed.jl | 1 + src/TransientDistributedCellField.jl | 2 +- 4 files changed, 151 insertions(+), 105 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index e0f93c99..10be3bc2 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -1,32 +1,37 @@ -abstract type DistributedCellDatum <: GridapType end - # DistributedCellPoint """ """ -struct DistributedCellPoint{A<:AbstractArray{<:CellPoint}} <: DistributedCellDatum +struct DistributedCellPoint{A<:AbstractArray{<:CellPoint},B<:DistributedTriangulation} <: CellDatum points::A + trian ::B end local_views(a::DistributedCellPoint) = a.points +CellData.get_triangulation(a::DistributedCellPoint) = a.trian # DistributedCellField """ """ -struct DistributedCellField{A,B} <: DistributedCellDatum +struct DistributedCellField{A,B,C} <: CellField fields::A - metadata::B + trian ::B + metadata::C function DistributedCellField( fields::AbstractArray{<:CellField}, + trian ::DistributedTriangulation, metadata=nothing) A = typeof(fields) - B = typeof(metadata) - new{A,B}(fields,metadata) + B = typeof(trian) + C = typeof(metadata) + new{A,B,C}(fields,trian,metadata) end end local_views(a::DistributedCellField) = a.fields +CellData.get_triangulation(a::DistributedCellField) = a.trian +CellData.DomainStyle(a::DistributedCellField) = DomainStyle(getany(a.fields)) # Constructors @@ -34,14 +39,14 @@ function CellData.CellField(f::Function,trian::DistributedTriangulation) fields = map(trian.trians) do t CellField(f,t) end - DistributedCellField(fields) + DistributedCellField(fields,trian) end function CellData.CellField(f::Number,trian::DistributedTriangulation) fields = map(trian.trians) do t CellField(f,t) end - DistributedCellField(fields) + DistributedCellField(fields,trian) end function CellData.CellField( @@ -49,7 +54,7 @@ function CellData.CellField( fields = map(f,trian.trians) do f,t CellField(f,t) end - DistributedCellField(fields) + DistributedCellField(fields,trian) end # Evaluation @@ -59,135 +64,168 @@ function (f::DistributedCellField)(x::DistributedCellPoint) end function Arrays.evaluate!(cache,f::DistributedCellField,x::DistributedCellPoint) - map(f.fields,x.points) do f,x + map(local_views(f),local_views(x)) do f,x evaluate!(nothing,f,x) end end # Operations +function _select_triangulation(fields,parents::DistributedCellField...) + trian_candidates = unique(objectid,map(get_triangulation,parents)) + _select_triangulation(fields,trian_candidates...) +end + +function _select_triangulation(fields,trian_candidates::DistributedTriangulation...) + if length(trian_candidates) == 1 + return first(trian_candidates) + end + + # Check if we can select one of the original triangulations + trians = map(local_views,trian_candidates) + t_id = map(fields,trians...) do f, trians + f_id = objectid(get_triangulation(f)) + return findfirst(tt -> objectid(tt) == f_id, trians) + end |> getany + if !isnothing(t_id) + return trians[t_id] + end + + # If not, check if we can build a new DistributedTriangulation based on one of the original models. + m_id = map(fields,trians...) do f, trians + f_id = objectid(get_background_model(get_triangulation(f))) + return findfirst(tt -> objectid(get_background_model(tt)) == f_id, trians) + end |> getany + return DistributedTriangulation(map(get_triangulation,fields),get_background_model(trians[m_id])) + + @error "Cannot select a triangulation for the operation" +end + function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField) fields = map(a.fields) do f evaluate!(nothing,k,f) end - DistributedCellField(fields) + DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!( cache,k::Operation,a::DistributedCellField,b::DistributedCellField) - fields = map(a.fields,b.fields) do f,g + fields = map(local_views(a),local_views(b)) do f,g evaluate!(nothing,k,f,g) end - DistributedCellField(fields) + trian = _select_triangulation(fields,a,b) + DistributedCellField(fields,trian) end function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::Number) fields = map(a.fields) do f evaluate!(nothing,k,f,b) end - DistributedCellField(fields) + DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!(cache,k::Operation,b::Number,a::DistributedCellField) fields = map(a.fields) do f evaluate!(nothing,k,b,f) end - DistributedCellField(fields) + DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::Function) fields = map(a.fields) do f evaluate!(nothing,k,f,b) end - DistributedCellField(fields) + DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!(cache,k::Operation,b::Function,a::DistributedCellField) fields = map(a.fields) do f evaluate!(nothing,k,b,f) end - DistributedCellField(fields) + DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField...) - fields = map(map(i->i.fields,a)...) do f... + fields = map(map(local_views,a)...) do f... evaluate!(nothing,k,f...) end - DistributedCellField(fields) + trian = _select_triangulation(fields,a...) + DistributedCellField(fields,trian) end # Composition -Base.:(∘)(f::Function,g::DistributedCellField) = Operation(f)(g) -Base.:(∘)(f::Function,g::Tuple{DistributedCellField,DistributedCellField}) = Operation(f)(g[1],g[2]) -Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Number}) = Operation(f)(g[1],g[2]) -Base.:(∘)(f::Function,g::Tuple{Number,DistributedCellField}) = Operation(f)(g[1],g[2]) -Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Function}) = Operation(f)(g[1],g[2]) -Base.:(∘)(f::Function,g::Tuple{Function,DistributedCellField}) = Operation(f)(g[1],g[2]) -Base.:(∘)(f::Function,g::Tuple{Vararg{DistributedCellField}}) = Operation(f)(g...) +# Base.:(∘)(f::Function,g::DistributedCellField) = Operation(f)(g) +# Base.:(∘)(f::Function,g::Tuple{DistributedCellField,DistributedCellField}) = Operation(f)(g[1],g[2]) +# Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Number}) = Operation(f)(g[1],g[2]) +# Base.:(∘)(f::Function,g::Tuple{Number,DistributedCellField}) = Operation(f)(g[1],g[2]) +# Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Function}) = Operation(f)(g[1],g[2]) +# Base.:(∘)(f::Function,g::Tuple{Function,DistributedCellField}) = Operation(f)(g[1],g[2]) +# Base.:(∘)(f::Function,g::Tuple{Vararg{DistributedCellField}}) = Operation(f)(g...) # Define some of the well known arithmetic ops # Unary ops -for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj) - @eval begin - ($op)(a::DistributedCellField) = Operation($op)(a) - end -end +#for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj) +# @eval begin +# ($op)(a::DistributedCellField) = Operation($op)(a) +# end +#end # Binary ops - -for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/) - @eval begin - ($op)(a::DistributedCellField,b::DistributedCellField) = Operation($op)(a,b) - ($op)(a::DistributedCellField,b::Number) = Operation($op)(a,b) - ($op)(a::Number,b::DistributedCellField) = Operation($op)(a,b) - ($op)(a::DistributedCellField,b::Function) = Operation($op)(a,b) - ($op)(a::Function,b::DistributedCellField) = Operation($op)(a,b) - end -end - -Base.broadcasted(f,a::DistributedCellField,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) -Base.broadcasted(f,a::Number,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) -Base.broadcasted(f,a::DistributedCellField,b::Number) = Operation((i,j)->f.(i,j))(a,b) -Base.broadcasted(f,a::Function,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) -Base.broadcasted(f,a::DistributedCellField,b::Function) = Operation((i,j)->f.(i,j))(a,b) -Base.broadcasted(::typeof(*),::typeof(∇),f::DistributedCellField) = Operation(Fields._extract_grad_diag)(∇(f)) -Base.broadcasted(::typeof(*),s::Fields.ShiftedNabla,f::DistributedCellField) = Operation(Fields._extract_grad_diag)(s(f)) - -dot(::typeof(∇),f::DistributedCellField) = divergence(f) -outer(::typeof(∇),f::DistributedCellField) = gradient(f) -outer(f::DistributedCellField,::typeof(∇)) = transpose(gradient(f)) -cross(::typeof(∇),f::DistributedCellField) = curl(f) +# +#for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/) +# @eval begin +# ($op)(a::DistributedCellField,b::DistributedCellField) = Operation($op)(a,b) +# ($op)(a::DistributedCellField,b::Number) = Operation($op)(a,b) +# ($op)(a::Number,b::DistributedCellField) = Operation($op)(a,b) +# ($op)(a::DistributedCellField,b::Function) = Operation($op)(a,b) +# ($op)(a::Function,b::DistributedCellField) = Operation($op)(a,b) +# end +#end +# +#Base.broadcasted(f,a::DistributedCellField,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) +#Base.broadcasted(f,a::Number,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) +#Base.broadcasted(f,a::DistributedCellField,b::Number) = Operation((i,j)->f.(i,j))(a,b) +#Base.broadcasted(f,a::Function,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) +#Base.broadcasted(f,a::DistributedCellField,b::Function) = Operation((i,j)->f.(i,j))(a,b) +#Base.broadcasted(::typeof(*),::typeof(∇),f::DistributedCellField) = Operation(Fields._extract_grad_diag)(∇(f)) +#Base.broadcasted(::typeof(*),s::Fields.ShiftedNabla,f::DistributedCellField) = Operation(Fields._extract_grad_diag)(s(f)) +# +#dot(::typeof(∇),f::DistributedCellField) = divergence(f) +#outer(::typeof(∇),f::DistributedCellField) = gradient(f) +#outer(f::DistributedCellField,::typeof(∇)) = transpose(gradient(f)) +#cross(::typeof(∇),f::DistributedCellField) = curl(f) # Differential ops function Fields.gradient(a::DistributedCellField) - DistributedCellField(map(gradient,a.fields)) + DistributedCellField(map(gradient,a.fields),get_triangulation(a)) end function Fields.divergence(a::DistributedCellField) - DistributedCellField(map(divergence,a.fields)) + DistributedCellField(map(divergence,a.fields),get_triangulation(a)) end function Fields.DIV(a::DistributedCellField) - DistributedCellField(map(DIV,a.fields)) + DistributedCellField(map(DIV,a.fields),get_triangulation(a)) end function Fields.∇∇(a::DistributedCellField) - DistributedCellField(map(∇∇,a.fields)) + DistributedCellField(map(∇∇,a.fields),get_triangulation(a)) end function Fields.curl(a::DistributedCellField) - DistributedCellField(map(curl,a.fields)) + DistributedCellField(map(curl,a.fields),get_triangulation(a)) end # Integration related """ """ -struct DistributedMeasure{A<:AbstractArray{<:Measure}} <: GridapType +struct DistributedMeasure{A<:AbstractArray{<:Measure},B<:DistributedTriangulation} <: GridapType measures::A + trian::B end local_views(a::DistributedMeasure) = a.measures @@ -196,18 +234,18 @@ function CellData.Measure(t::DistributedTriangulation,args...) measures = map(t.trians) do trian Measure(trian,args...) end - DistributedMeasure(measures) + DistributedMeasure(measures,t) end function CellData.Measure(tt::DistributedTriangulation{Dc,Dp},it::DistributedTriangulation{Dc,Dp},args...) where {Dc,Dp} measures = map(local_views(tt),local_views(it)) do ttrian, itrian Measure(ttrian,itrian,args...) end - return DistributedMeasure(measures) + return DistributedMeasure(measures,it) end function CellData.get_cell_points(a::DistributedMeasure) - DistributedCellPoint(map(get_cell_points,a.measures)) + DistributedCellPoint(map(get_cell_points,a.measures),a.trian) end """ @@ -275,28 +313,30 @@ end # Triangulation related function CellData.get_cell_points(a::DistributedTriangulation) - DistributedCellPoint(map(get_cell_points,a.trians)) + DistributedCellPoint(map(get_cell_points,a.trians),a) end function CellData.get_normal_vector(a::DistributedTriangulation) fields = map(get_normal_vector,a.trians) - DistributedCellField(fields) + DistributedCellField(fields,a) end # Skeleton related -function DistributedCellField(a::AbstractArray{<:SkeletonPair}) +function DistributedCellField(a::AbstractArray{<:SkeletonPair},trian::DistributedTriangulation) plus, minus = map(s->(s.plus,s.minus),a) |> tuple_of_arrays - dplus = DistributedCellField(plus) - dminus = DistributedCellField(minus) + tplus = _select_triangulation(plus,trian) + tminus = _select_triangulation(minus,trian) + dplus = DistributedCellField(plus,tplus) + dminus = DistributedCellField(minus,tminus) SkeletonPair(dplus,dminus) end function Base.getproperty(x::DistributedCellField, sym::Symbol) if sym in (:⁺,:plus) - DistributedCellField(map(i->i.plus,x.fields)) + DistributedCellField(map(i->i.plus,local_views(x)),get_triangulation(x)) elseif sym in (:⁻, :minus) - DistributedCellField(map(i->i.minus,x.fields)) + DistributedCellField(map(i->i.minus,local_views(x)),get_triangulation(x)) else getfield(x, sym) end @@ -325,19 +365,16 @@ function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellFi SkeletonPair(plus,minus) end -CellData.jump(a::DistributedCellField) = DistributedCellField(map(jump,a.fields)) +CellData.jump(a::DistributedCellField) = DistributedCellField(map(jump,a.fields),get_triangulation(a)) CellData.jump(a::SkeletonPair{<:DistributedCellField}) = a.⁺ + a.⁻ -CellData.mean(a::DistributedCellField) = DistributedCellField(map(mean,a.fields)) +CellData.mean(a::DistributedCellField) = DistributedCellField(map(mean,a.fields),get_triangulation(a)) # DistributedCellDof -struct DistributedCellDof{A} <: DistributedCellDatum +struct DistributedCellDof{A,B} <: CellDatum dofs::A - function DistributedCellDof(dofs::AbstractArray{<:CellDof}) - A = typeof(dofs) - new{A}(dofs) - end + trian::B end local_views(s::DistributedCellDof) = s.dofs @@ -351,10 +388,10 @@ function Gridap.Arrays.evaluate!(cache,s::DistributedCellDof,f::DistributedCellF end function Gridap.Arrays.evaluate!(cache, ::DistributedCellField, ::DistributedCellDof) -@unreachable """\n -CellField (f) objects cannot be evaluated at CellDof (s) objects. -However, CellDofs objects can be evaluated at CellField objects. -Did you mean evaluate(f,s) instead of evaluate(s,f), i.e. -f(s) instead of s(f)? -""" + @unreachable """\n + CellField (f) objects cannot be evaluated at CellDof (s) objects. + However, CellDofs objects can be evaluated at CellField objects. + Did you mean evaluate(f,s) instead of evaluate(s,f), i.e. + f(s) instead of s(f)? + """ end \ No newline at end of file diff --git a/src/FESpaces.jl b/src/FESpaces.jl index cd922ddd..8bed308a 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -130,8 +130,6 @@ function generate_gids( cell_to_ldofs::AbstractArray{<:AbstractArray}, nldofs::AbstractArray{<:Integer}) - ngcells = length(cell_range) - # Find and count number owned dofs ldof_to_owner, nodofs = map(partition(cell_range),cell_to_ldofs,nldofs) do indices,cell_to_ldofs,nldofs ldof_to_owner = fill(Int32(0),nldofs) @@ -260,7 +258,7 @@ end """ """ -const DistributedSingleFieldFEFunction = DistributedCellField{A,<:DistributedFEFunctionData{T}} where {A,T} +const DistributedSingleFieldFEFunction = DistributedCellField{A,B,<:DistributedFEFunctionData{T}} where {A,B,T} function FESpaces.get_free_dof_values(uh::DistributedSingleFieldFEFunction) uh.metadata.free_values @@ -269,21 +267,25 @@ end # Single field related """ """ -struct DistributedSingleFieldFESpace{A,B,C} <: DistributedFESpace +struct DistributedSingleFieldFESpace{A,B,C,D} <: DistributedFESpace spaces::A gids::B - vector_type::Type{C} + trian::C + vector_type::Type{D} function DistributedSingleFieldFESpace( spaces::AbstractArray{<:SingleFieldFESpace}, gids::PRange, - vector_type::Type{C}) where C + trian::DistributedTriangulation, + vector_type::Type{D}) where D A = typeof(spaces) B = typeof(gids) - new{A,B,C}(spaces,gids,vector_type) + C = typeof(trian) + new{A,B,C,D}(spaces,gids,trian,vector_type) end end local_views(a::DistributedSingleFieldFESpace) = a.spaces +CellData.get_triangulation(a::DistributedSingleFieldFESpace) = a.trian function FESpaces.get_vector_type(fs::DistributedSingleFieldFESpace) fs.vector_type @@ -338,8 +340,9 @@ function _EvaluationFunction(func, f::DistributedSingleFieldFESpace,x::AbstractVector,isconsistent=false) free_values = change_ghost(x,f.gids,is_consistent=isconsistent,make_consistent=true) fields = map(func,f.spaces,partition(free_values)) + trian = get_triangulation(f) metadata = DistributedFEFunctionData(free_values) - DistributedCellField(fields,metadata) + DistributedCellField(fields,trian,metadata) end function _EvaluationFunction(func, @@ -347,56 +350,60 @@ function _EvaluationFunction(func, dirichlet_values::AbstractArray{<:AbstractVector},isconsistent=false) free_values = change_ghost(x,f.gids,is_consistent=isconsistent,make_consistent=true) fields = map(func,f.spaces,partition(free_values),dirichlet_values) + trian = get_triangulation(f) metadata = DistributedFEFunctionData(free_values) - DistributedCellField(fields,metadata) + DistributedCellField(fields,trian,metadata) end function FESpaces.get_fe_basis(f::DistributedSingleFieldFESpace) fields = map(get_fe_basis,f.spaces) - DistributedCellField(fields) + trian = get_triangulation(f) + DistributedCellField(fields,trian) end function FESpaces.get_trial_fe_basis(f::DistributedSingleFieldFESpace) fields = map(get_trial_fe_basis,f.spaces) - DistributedCellField(fields) + trian = get_triangulation(f) + DistributedCellField(fields,trian) end function FESpaces.get_fe_dof_basis(f::DistributedSingleFieldFESpace) - dofs = map(get_fe_dof_basis,local_views(f)) - DistributedCellDof(dofs) + dofs = map(get_fe_dof_basis,local_views(f)) + trian = get_triangulation(f) + DistributedCellDof(dofs,trian) end function FESpaces.TrialFESpace(f::DistributedSingleFieldFESpace) spaces = map(TrialFESpace,f.spaces) - DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type) + DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type) end function FESpaces.TrialFESpace(f::DistributedSingleFieldFESpace,fun) spaces = map(f.spaces) do s TrialFESpace(s,fun) end - DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type) + DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type) end function FESpaces.TrialFESpace(fun,f::DistributedSingleFieldFESpace) spaces = map(f.spaces) do s TrialFESpace(fun,s) end - DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type) + DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type) end function FESpaces.TrialFESpace!(f::DistributedSingleFieldFESpace,fun) spaces = map(f.spaces) do s TrialFESpace!(s,fun) end - DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type) + DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type) end function FESpaces.HomogeneousTrialFESpace(f::DistributedSingleFieldFESpace) spaces = map(f.spaces) do s HomogeneousTrialFESpace(s) end - DistributedSingleFieldFESpace(spaces,f.gids,f.vector_type) + DistributedSingleFieldFESpace(spaces,f.gids,f.trian,f.vector_type) end function generate_gids( @@ -478,8 +485,9 @@ function FESpaces.FESpace(model::DistributedDiscreteModel,reffe;split_own_and_gh FESpace(m,reffe;kwargs...) end gids = generate_gids(model,spaces) + trian = DistributedTriangulation(map(get_triangulation,spaces),model) vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost) - DistributedSingleFieldFESpace(spaces,gids,vector_type) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) end function FESpaces.FESpace(_trian::DistributedTriangulation,reffe;split_own_and_ghost=false,kwargs...) @@ -492,7 +500,7 @@ function FESpaces.FESpace(_trian::DistributedTriangulation,reffe;split_own_and_g nldofs = map(num_free_dofs,spaces) gids = generate_gids(trian_gids,cell_to_ldofs,nldofs) vector_type = _find_vector_type(spaces,gids;split_own_and_ghost=split_own_and_ghost) - DistributedSingleFieldFESpace(spaces,gids,vector_type) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) end function _find_vector_type(spaces,gids;split_own_and_ghost=false) diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl index 93e1077c..45560fe7 100644 --- a/src/GridapDistributed.jl +++ b/src/GridapDistributed.jl @@ -17,6 +17,7 @@ using Gridap.ODEs using MPI using PartitionedArrays const PArrays = PartitionedArrays +using PartitionedArrays: getany using SparseArrays using WriteVTK diff --git a/src/TransientDistributedCellField.jl b/src/TransientDistributedCellField.jl index 3186d986..b8a63709 100644 --- a/src/TransientDistributedCellField.jl +++ b/src/TransientDistributedCellField.jl @@ -1,5 +1,5 @@ # Transient Distributed CellField -abstract type TransientDistributedCellField <: DistributedCellDatum end +abstract type TransientDistributedCellField <: CellDatum end # Transient SingleField struct TransientSingleFieldDistributedCellField{A} <: TransientDistributedCellField From e7651c163f24bee99c1a6ce58b7f98c0abb2772d Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 19 Dec 2023 14:38:59 +1100 Subject: [PATCH 03/25] Some bugfixes --- src/CellData.jl | 11 +++++++---- src/DivConformingFESpaces.jl | 3 ++- src/GridapDistributed.jl | 6 +++--- 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index 10be3bc2..cc05fcf5 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -83,20 +83,23 @@ function _select_triangulation(fields,trian_candidates::DistributedTriangulation # Check if we can select one of the original triangulations trians = map(local_views,trian_candidates) - t_id = map(fields,trians...) do f, trians + t_id = map(fields,trians...) do f, trians... f_id = objectid(get_triangulation(f)) return findfirst(tt -> objectid(tt) == f_id, trians) end |> getany if !isnothing(t_id) - return trians[t_id] + return trian_candidates[t_id] end # If not, check if we can build a new DistributedTriangulation based on one of the original models. - m_id = map(fields,trians...) do f, trians + m_id = map(fields,trians...) do f, trians... f_id = objectid(get_background_model(get_triangulation(f))) return findfirst(tt -> objectid(get_background_model(tt)) == f_id, trians) end |> getany - return DistributedTriangulation(map(get_triangulation,fields),get_background_model(trians[m_id])) + if !isnothing(m_id) + model = get_background_model(trian_candidates[m_id]) + return DistributedTriangulation(map(get_triangulation,fields),model) + end @error "Cannot select a triangulation for the operation" end diff --git a/src/DivConformingFESpaces.jl b/src/DivConformingFESpaces.jl index ac31e505..c814e75e 100644 --- a/src/DivConformingFESpaces.jl +++ b/src/DivConformingFESpaces.jl @@ -27,8 +27,9 @@ function _common_fe_space_constructor(model,cell_reffes;conformity,kwargs...) FESpace(m, cell_fe; kwargs...) end gids = generate_gids(model,spaces) + trian = DistributedTriangulation(map(get_triangulation,spaces),model) vector_type = _find_vector_type(spaces,gids) - DistributedSingleFieldFESpace(spaces,gids,vector_type) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) end function _generate_sign_flips(model,cell_reffes) diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl index 45560fe7..f01581de 100644 --- a/src/GridapDistributed.jl +++ b/src/GridapDistributed.jl @@ -57,11 +57,11 @@ include("DivConformingFESpaces.jl") include("MultiField.jl") -include("TransientDistributedCellField.jl") +#include("TransientDistributedCellField.jl") -include("TransientMultiFieldDistributedCellField.jl") +#include("TransientMultiFieldDistributedCellField.jl") -include("TransientFESpaces.jl") +#include("TransientFESpaces.jl") include("Adaptivity.jl") From 3c33c5d951d6a3c33208ec2e6938cc66062f61c6 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 19 Dec 2023 14:54:53 +1100 Subject: [PATCH 04/25] Removed unnecessary methods --- src/CellData.jl | 92 +++++++------------------------------------------ 1 file changed, 13 insertions(+), 79 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index cc05fcf5..2b7eda37 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -9,6 +9,7 @@ end local_views(a::DistributedCellPoint) = a.points CellData.get_triangulation(a::DistributedCellPoint) = a.trian +CellData.DomainStyle(a::DistributedCellPoint) = DomainStyle(getany(local_views(a))) # DistributedCellField """ @@ -31,7 +32,7 @@ end local_views(a::DistributedCellField) = a.fields CellData.get_triangulation(a::DistributedCellField) = a.trian -CellData.DomainStyle(a::DistributedCellField) = DomainStyle(getany(a.fields)) +CellData.DomainStyle(a::DistributedCellField) = DomainStyle(getany(local_views(a))) # Constructors @@ -69,8 +70,8 @@ function Arrays.evaluate!(cache,f::DistributedCellField,x::DistributedCellPoint) end end -# Operations - +# Given local CellFields and a set of original DistributedTriangulations, +# returns the DistributedTriangulation where the local CellFields are defined. function _select_triangulation(fields,parents::DistributedCellField...) trian_candidates = unique(objectid,map(get_triangulation,parents)) _select_triangulation(fields,trian_candidates...) @@ -104,8 +105,10 @@ function _select_triangulation(fields,trian_candidates::DistributedTriangulation @error "Cannot select a triangulation for the operation" end +# Operations + function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField) - fields = map(a.fields) do f + fields = map(local_views(a)) do f evaluate!(nothing,k,f) end DistributedCellField(fields,get_triangulation(a)) @@ -121,28 +124,28 @@ function Arrays.evaluate!( end function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::Number) - fields = map(a.fields) do f + fields = map(local_views(a)) do f evaluate!(nothing,k,f,b) end DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!(cache,k::Operation,b::Number,a::DistributedCellField) - fields = map(a.fields) do f + fields = map(local_views(a)) do f evaluate!(nothing,k,b,f) end DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::Function) - fields = map(a.fields) do f + fields = map(local_views(a)) do f evaluate!(nothing,k,f,b) end DistributedCellField(fields,get_triangulation(a)) end function Arrays.evaluate!(cache,k::Operation,b::Function,a::DistributedCellField) - fields = map(a.fields) do f + fields = map(local_views(a)) do f evaluate!(nothing,k,b,f) end DistributedCellField(fields,get_triangulation(a)) @@ -156,51 +159,6 @@ function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField...) DistributedCellField(fields,trian) end -# Composition - -# Base.:(∘)(f::Function,g::DistributedCellField) = Operation(f)(g) -# Base.:(∘)(f::Function,g::Tuple{DistributedCellField,DistributedCellField}) = Operation(f)(g[1],g[2]) -# Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Number}) = Operation(f)(g[1],g[2]) -# Base.:(∘)(f::Function,g::Tuple{Number,DistributedCellField}) = Operation(f)(g[1],g[2]) -# Base.:(∘)(f::Function,g::Tuple{DistributedCellField,Function}) = Operation(f)(g[1],g[2]) -# Base.:(∘)(f::Function,g::Tuple{Function,DistributedCellField}) = Operation(f)(g[1],g[2]) -# Base.:(∘)(f::Function,g::Tuple{Vararg{DistributedCellField}}) = Operation(f)(g...) - -# Define some of the well known arithmetic ops - -# Unary ops - -#for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj) -# @eval begin -# ($op)(a::DistributedCellField) = Operation($op)(a) -# end -#end - -# Binary ops -# -#for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/) -# @eval begin -# ($op)(a::DistributedCellField,b::DistributedCellField) = Operation($op)(a,b) -# ($op)(a::DistributedCellField,b::Number) = Operation($op)(a,b) -# ($op)(a::Number,b::DistributedCellField) = Operation($op)(a,b) -# ($op)(a::DistributedCellField,b::Function) = Operation($op)(a,b) -# ($op)(a::Function,b::DistributedCellField) = Operation($op)(a,b) -# end -#end -# -#Base.broadcasted(f,a::DistributedCellField,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) -#Base.broadcasted(f,a::Number,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) -#Base.broadcasted(f,a::DistributedCellField,b::Number) = Operation((i,j)->f.(i,j))(a,b) -#Base.broadcasted(f,a::Function,b::DistributedCellField) = Operation((i,j)->f.(i,j))(a,b) -#Base.broadcasted(f,a::DistributedCellField,b::Function) = Operation((i,j)->f.(i,j))(a,b) -#Base.broadcasted(::typeof(*),::typeof(∇),f::DistributedCellField) = Operation(Fields._extract_grad_diag)(∇(f)) -#Base.broadcasted(::typeof(*),s::Fields.ShiftedNabla,f::DistributedCellField) = Operation(Fields._extract_grad_diag)(s(f)) -# -#dot(::typeof(∇),f::DistributedCellField) = divergence(f) -#outer(::typeof(∇),f::DistributedCellField) = gradient(f) -#outer(f::DistributedCellField,::typeof(∇)) = transpose(gradient(f)) -#cross(::typeof(∇),f::DistributedCellField) = curl(f) - # Differential ops function Fields.gradient(a::DistributedCellField) @@ -345,34 +303,9 @@ function Base.getproperty(x::DistributedCellField, sym::Symbol) end end -function Base.propertynames(x::DistributedCellField, private::Bool=false) - (fieldnames(typeof(x))...,:⁺,:plus,:⁻,:minus) -end - -for op in (:outer,:*,:dot) - @eval begin - ($op)(a::DistributedCellField,b::SkeletonPair{<:DistributedCellField}) = Operation($op)(a,b) - ($op)(a::SkeletonPair{<:DistributedCellField},b::DistributedCellField) = Operation($op)(a,b) - end -end - -function Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::SkeletonPair{<:DistributedCellField}) - plus = k(a.plus,b.plus) - minus = k(a.minus,b.minus) - SkeletonPair(plus,minus) -end - -function Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellField},b::DistributedCellField) - plus = k(a.plus,b.plus) - minus = k(a.minus,b.minus) - SkeletonPair(plus,minus) -end - CellData.jump(a::DistributedCellField) = DistributedCellField(map(jump,a.fields),get_triangulation(a)) -CellData.jump(a::SkeletonPair{<:DistributedCellField}) = a.⁺ + a.⁻ CellData.mean(a::DistributedCellField) = DistributedCellField(map(mean,a.fields),get_triangulation(a)) - # DistributedCellDof struct DistributedCellDof{A,B} <: CellDatum @@ -381,12 +314,13 @@ struct DistributedCellDof{A,B} <: CellDatum end local_views(s::DistributedCellDof) = s.dofs +CellData.get_triangulation(s::DistributedCellDof) = s.trian (a::DistributedCellDof)(f) = evaluate(a,f) function Gridap.Arrays.evaluate!(cache,s::DistributedCellDof,f::DistributedCellField) map(local_views(s),local_views(f)) do s, f - evaluate!(nothing,s,f) + evaluate!(nothing,s,f) end end From 79f2dc8fafb68c95f6fe10b44bf3bbbe4aa3e6d8 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 19 Dec 2023 18:10:07 +1100 Subject: [PATCH 05/25] DMultiFieldFEFunction -> DMultiFieldCellField Now the distributed multifield implementation follows a bit more what the single-field is doing, and brings us closer to Gridap. I have renamed DMultifieldFEFunction as DMultiFieldCellField, and have created aliases for the FEFunction and FEBasis specialisations, --- src/CellData.jl | 16 +++- src/FESpaces.jl | 2 +- src/GridapDistributed.jl | 6 +- src/MultiField.jl | 131 +++++++++++++++------------ src/ODEs.jl | 37 ++++++++ src/TransientDistributedCellField.jl | 107 ---------------------- src/TransientFESpaces.jl | 44 --------- 7 files changed, 124 insertions(+), 219 deletions(-) create mode 100644 src/ODEs.jl delete mode 100644 src/TransientDistributedCellField.jl diff --git a/src/CellData.jl b/src/CellData.jl index 2b7eda37..68983440 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -9,7 +9,10 @@ end local_views(a::DistributedCellPoint) = a.points CellData.get_triangulation(a::DistributedCellPoint) = a.trian -CellData.DomainStyle(a::DistributedCellPoint) = DomainStyle(getany(local_views(a))) + +function CellData.DomainStyle(::Type{<:DistributedCellPoint{A}}) where A + DomainStyle(eltype(A)) +end # DistributedCellField """ @@ -32,7 +35,10 @@ end local_views(a::DistributedCellField) = a.fields CellData.get_triangulation(a::DistributedCellField) = a.trian -CellData.DomainStyle(a::DistributedCellField) = DomainStyle(getany(local_views(a))) + +function CellData.DomainStyle(::Type{<:DistributedCellField{A}}) where A + DomainStyle(eltype(A)) +end # Constructors @@ -308,7 +314,7 @@ CellData.mean(a::DistributedCellField) = DistributedCellField(map(mean,a.fields) # DistributedCellDof -struct DistributedCellDof{A,B} <: CellDatum +struct DistributedCellDof{A<:AbstractArray{<:CellDof},B<:DistributedTriangulation} <: CellDatum dofs::A trian::B end @@ -316,6 +322,10 @@ end local_views(s::DistributedCellDof) = s.dofs CellData.get_triangulation(s::DistributedCellDof) = s.trian +function CellData.DomainStyle(::Type{<:DistributedCellDof{A}}) where A + DomainStyle(eltype(A)) +end + (a::DistributedCellDof)(f) = evaluate(a,f) function Gridap.Arrays.evaluate!(cache,s::DistributedCellDof,f::DistributedCellField) diff --git a/src/FESpaces.jl b/src/FESpaces.jl index 8bed308a..8788d099 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -258,7 +258,7 @@ end """ """ -const DistributedSingleFieldFEFunction = DistributedCellField{A,B,<:DistributedFEFunctionData{T}} where {A,B,T} +const DistributedSingleFieldFEFunction{A,B,T} = DistributedCellField{A,B,DistributedFEFunctionData{T}} function FESpaces.get_free_dof_values(uh::DistributedSingleFieldFEFunction) uh.metadata.free_values diff --git a/src/GridapDistributed.jl b/src/GridapDistributed.jl index f01581de..0472ef61 100644 --- a/src/GridapDistributed.jl +++ b/src/GridapDistributed.jl @@ -57,11 +57,7 @@ include("DivConformingFESpaces.jl") include("MultiField.jl") -#include("TransientDistributedCellField.jl") - -#include("TransientMultiFieldDistributedCellField.jl") - -#include("TransientFESpaces.jl") +include("ODEs.jl") include("Adaptivity.jl") diff --git a/src/MultiField.jl b/src/MultiField.jl index 2d535ab4..40f61d51 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -1,31 +1,55 @@ -struct DistributedMultiFieldFEFunction{A,B,C} <: GridapType +# DistributedMultiFieldCellField +struct DistributedMultiFieldCellField{A,B,C} <: CellField field_fe_fun::A part_fe_fun::B - free_values::C - function DistributedMultiFieldFEFunction( - field_fe_fun::AbstractVector{<:DistributedSingleFieldFEFunction}, - part_fe_fun::AbstractArray{<:MultiFieldFEFunction}, - free_values::AbstractVector) + metadata::C + function DistributedMultiFieldCellField( + field_fe_fun::AbstractVector{<:DistributedCellField}, + part_fe_fun::AbstractArray{<:CellField}, + metadata=nothing) A = typeof(field_fe_fun) B = typeof(part_fe_fun) - C = typeof(free_values) - new{A,B,C}(field_fe_fun,part_fe_fun,free_values) + C = typeof(metadata) + new{A,B,C}(field_fe_fun,part_fe_fun,metadata) end end +function CellData.get_triangulation(a::DistributedMultiFieldCellField) + trians = map(get_triangulation,a.field_fe_fun) + @check all(map(t -> t === first(trians), trians)) + return first(trians) +end -function FESpaces.get_free_dof_values(uh::DistributedMultiFieldFEFunction) - uh.free_values +function CellData.DomainStyle(::Type{<:DistributedMultiFieldCellField{A,B}}) where {A,B} + DomainStyle(eltype(B)) +end + +local_views(a::DistributedMultiFieldCellField) = a.part_fe_fun +local_views(a::Vector{<:DistributedCellField}) = map(local_views,a) + +MultiField.num_fields(m::DistributedMultiFieldCellField) = length(m.field_fe_fun) +Base.iterate(m::DistributedMultiFieldCellField) = iterate(m.field_fe_fun) +Base.iterate(m::DistributedMultiFieldCellField,state) = iterate(m.field_fe_fun,state) +Base.getindex(m::DistributedMultiFieldCellField,field_id::Integer) = m.field_fe_fun[field_id] + +# DistributedMultiFieldFEFunction + +const DistributedMultiFieldFEFunction{A,B,T} = DistributedMultiFieldCellField{A,B,DistributedFEFunctionData{T}} + +function DistributedMultiFieldFEFunction( + field_fe_fun::AbstractVector{<:DistributedSingleFieldFEFunction}, + part_fe_fun::AbstractArray{<:MultiFieldFEFunction}, + free_values::AbstractVector) + metadata = DistributedFEFunctionData(free_values) + DistributedMultiFieldCellField(field_fe_fun,part_fe_fun,metadata) end -local_views(a::DistributedMultiFieldFEFunction) = a.part_fe_fun -MultiField.num_fields(m::DistributedMultiFieldFEFunction) = length(m.field_fe_fun) -Base.iterate(m::DistributedMultiFieldFEFunction) = iterate(m.field_fe_fun) -Base.iterate(m::DistributedMultiFieldFEFunction,state) = iterate(m.field_fe_fun,state) -Base.getindex(m::DistributedMultiFieldFEFunction,field_id::Integer) = m.field_fe_fun[field_id] +function FESpaces.get_free_dof_values(uh::DistributedMultiFieldFEFunction) + uh.metadata.free_values +end -local_views(a::Vector{<:DistributedCellField}) = [ai.fields for ai in a] +# DistributedMultiFieldFESpace """ """ @@ -47,6 +71,14 @@ struct DistributedMultiFieldFESpace{MS,A,B,C,D} <: DistributedFESpace end end +function CellData.get_triangulation(a::DistributedMultiFieldFESpace) + trians = map(get_triangulation,a.field_fe_space) + @check all(map(t -> t === first(trians), trians)) + return first(trians) +end + +MultiField.MultiFieldStyle(::Type{<:DistributedMultiFieldFESpace{MS}}) where MS = MS() + local_views(a::DistributedMultiFieldFESpace) = a.part_fe_space MultiField.num_fields(m::DistributedMultiFieldFESpace) = length(m.field_fe_space) Base.iterate(m::DistributedMultiFieldFESpace) = iterate(m.field_fe_space) @@ -178,49 +210,6 @@ function FESpaces.interpolate_everywhere( DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) end - -""" -""" -struct DistributedMultiFieldFEBasis{A,B} <: GridapType - field_fe_basis::A - part_fe_basis::B - function DistributedMultiFieldFEBasis( - field_fe_basis::AbstractVector{<:DistributedCellField}, - part_fe_basis::AbstractArray{<:MultiFieldCellField}) - A = typeof(field_fe_basis) - B = typeof(part_fe_basis) - new{A,B}(field_fe_basis,part_fe_basis) - end -end - -local_views(a::DistributedMultiFieldFEBasis) = a.part_fe_basis -MultiField.num_fields(m::DistributedMultiFieldFEBasis) = length(m.field_fe_basis) -Base.iterate(m::DistributedMultiFieldFEBasis) = iterate(m.field_fe_basis) -Base.iterate(m::DistributedMultiFieldFEBasis,state) = iterate(m.field_fe_basis,state) -Base.getindex(m::DistributedMultiFieldFEBasis,field_id::Integer) = m.field_fe_basis[field_id] - -function FESpaces.get_fe_basis(f::DistributedMultiFieldFESpace) - part_mbasis = map(get_fe_basis,f.part_fe_space) - field_fe_basis = DistributedCellField[] - for i in 1:num_fields(f) - basis_i = map(b->b[i],part_mbasis) - bi = DistributedCellField(basis_i) - push!(field_fe_basis,bi) - end - DistributedMultiFieldFEBasis(field_fe_basis,part_mbasis) -end - -function FESpaces.get_trial_fe_basis(f::DistributedMultiFieldFESpace) - part_mbasis = map(get_trial_fe_basis,f.part_fe_space) - field_fe_basis = DistributedCellField[] - for i in 1:num_fields(f) - basis_i = map(b->b[i],part_mbasis) - bi = DistributedCellField(basis_i) - push!(field_fe_basis,bi) - end - DistributedMultiFieldFEBasis(field_fe_basis,part_mbasis) -end - function FESpaces.TrialFESpace(objects,a::DistributedMultiFieldFESpace) TrialFESpace(a,objects) end @@ -237,6 +226,30 @@ function FESpaces.TrialFESpace(a::DistributedMultiFieldFESpace{MS},objects) wher DistributedMultiFieldFESpace(f_dspace,p_mspace,gids,vector_type) end +# DistributedMultiFieldFEBasis + +const DistributedMultiFieldFEBasis{A,B<:AbstractArray{<:FEBasis}} = DistributedMultiFieldCellField{A,B} + +function FESpaces.get_fe_basis(f::DistributedMultiFieldFESpace) + part_mbasis = map(get_fe_basis,f.part_fe_space) + field_fe_basis = map(1:num_fields(f)) do i + space_i = f.field_fe_space[i] + basis_i = map(b->b[i],part_mbasis) + DistributedCellField(basis_i,get_triangulation(space_i)) + end + DistributedMultiFieldCellField(field_fe_basis,part_mbasis) +end + +function FESpaces.get_trial_fe_basis(f::DistributedMultiFieldFESpace) + part_mbasis = map(get_trial_fe_basis,f.part_fe_space) + field_fe_basis = map(1:num_fields(f)) do i + space_i = f.field_fe_space[i] + basis_i = map(b->b[i],part_mbasis) + DistributedCellField(basis_i,get_triangulation(space_i)) + end + DistributedMultiFieldCellField(field_fe_basis,part_mbasis) +end + # Factory function MultiField.MultiFieldFESpace( diff --git a/src/ODEs.jl b/src/ODEs.jl new file mode 100644 index 00000000..d94ddd80 --- /dev/null +++ b/src/ODEs.jl @@ -0,0 +1,37 @@ + +# SingleField + +function ODEs.TransientCellField(f::DistributedCellField,derivatives::Tuple) + fields = map(local_views(f),map(local_views,derivatives)...) do f, derivatives... + ODEs.TransientCellField(f,Tuple(derivatives)) + end + DistributedCellField(fields,get_triangulation(f),f.metadata) +end + +function ODEs.time_derivative(f::DistributedCellField) + fields = map(local_views(f)) do field + ODEs.time_derivative(field) + end + DistributedCellField(fields,get_triangulation(f)) +end + +# MultiField + +function ODEs.TransientCellField(f::DistributedMultiFieldFEFunction,derivatives::Tuple) + field_fe_fun = map(1:num_fields(f)) do i + f_i = f[i] + df_i = Tuple(map(df -> df[i],derivatives)) + ODEs.TransientCellField(f_i,df_i) + end + part_fe_fun = map(local_views(f),map(local_views,derivatives)...) do f, derivatives... + ODEs.TransientCellField(f,Tuple(derivatives)) + end + free_values = get_free_dof_values(f) + DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) +end + +function ODEs.time_derivative(f::DistributedMultiFieldCellField) + field_fe_fun = map(ODEs.time_derivative,f.field_fe_fun) + part_fe_fun = map(ODEs.time_derivative,f.part_fe_fun) + DistributedMultiFieldCellField(field_fe_fun,part_fe_fun) +end diff --git a/src/TransientDistributedCellField.jl b/src/TransientDistributedCellField.jl deleted file mode 100644 index b8a63709..00000000 --- a/src/TransientDistributedCellField.jl +++ /dev/null @@ -1,107 +0,0 @@ -# Transient Distributed CellField -abstract type TransientDistributedCellField <: CellDatum end - -# Transient SingleField -struct TransientSingleFieldDistributedCellField{A} <: TransientDistributedCellField - cellfield::A - derivatives::Tuple -end - -local_views(f::TransientSingleFieldDistributedCellField) = local_views(f.cellfield) - -# Constructors -function ODEs.TransientCellField(single_field::DistributedSingleFieldFEFunction,derivatives::Tuple) - TransientSingleFieldDistributedCellField(single_field,derivatives) -end - -function ODEs.TransientCellField(single_field::DistributedCellField,derivatives::Tuple) - TransientSingleFieldDistributedCellField(single_field,derivatives) -end - -# Time derivative -function ODEs.∂t(f::TransientDistributedCellField) - cellfield, derivatives = first_and_tail(f.derivatives) - TransientCellField(cellfield,derivatives) -end - -ODEs.∂tt(f::TransientDistributedCellField) = ∂t(∂t(f)) - -# Integration related -function Fields.integrate(f::TransientDistributedCellField,b::DistributedMeasure) - integrate(f.cellfield,b) -end - -# Differential Operations -Fields.gradient(f::TransientDistributedCellField) = gradient(f.cellfield) -Fields.∇∇(f::TransientDistributedCellField) = ∇∇(f.cellfield) - -# Unary ops -for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj) - @eval begin - ($op)(a::TransientDistributedCellField) = ($op)(a.cellfield) - end -end - -# Binary ops -for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/) - @eval begin - ($op)(a::TransientDistributedCellField,b::TransientDistributedCellField) = ($op)(a.cellfield,b.cellfield) - ($op)(a::TransientDistributedCellField,b::DistributedCellField) = ($op)(a.cellfield,b) - ($op)(a::DistributedCellField,b::TransientDistributedCellField) = ($op)(a,b.cellfield) - ($op)(a::TransientDistributedCellField,b::Number) = ($op)(a.cellfield,b) - ($op)(a::Number,b::TransientDistributedCellField) = ($op)(a,b.cellfield) - ($op)(a::TransientDistributedCellField,b::Function) = ($op)(a.cellfield,b) - ($op)(a::Function,b::TransientDistributedCellField) = ($op)(a,b.cellfield) - end -end - -Base.broadcasted(f,a::TransientDistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a.cellfield,b.cellfield) -Base.broadcasted(f,a::TransientDistributedCellField,b::DistributedCellField) = broadcasted(f,a.cellfield,b) -Base.broadcasted(f,a::DistributedCellField,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield) -Base.broadcasted(f,a::Number,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield) -Base.broadcasted(f,a::TransientDistributedCellField,b::Number) = broadcasted(f,a.cellfield,b) -Base.broadcasted(f,a::Function,b::TransientDistributedCellField) = broadcasted(f,a,b.cellfield) -Base.broadcasted(f,a::TransientDistributedCellField,b::Function) = broadcasted(f,a.cellfield,b) -Base.broadcasted(a::typeof(*),b::typeof(∇),f::TransientDistributedCellField) = broadcasted(a,b,f.cellfield) -Base.broadcasted(a::typeof(*),s::Fields.ShiftedNabla,f::TransientDistributedCellField) = broadcasted(a,s,f.cellfield) - -dot(::typeof(∇),f::TransientDistributedCellField) = dot(∇,f.cellfield) -outer(::typeof(∇),f::TransientDistributedCellField) = outer(∇,f.cellfield) -outer(f::TransientDistributedCellField,::typeof(∇)) = outer(f.cellfield,∇) -cross(::typeof(∇),f::TransientDistributedCellField) = cross(∇,f.cellfield) - -Gridap.Arrays.evaluate!(cache,k::Operation,a::TransientDistributedCellField,b::DistributedCellField) = evaluate!(cache,k,a.cellfield,b) -Gridap.Arrays.evaluate!(cache,k::Operation,a::DistributedCellField,b::TransientDistributedCellField) = evaluate!(cache,k,a,b.cellfield) - -Base.:(∘)(f::Function,g::Tuple{TransientDistributedCellField,DistributedCellField}) = Operation(f)(g[1],g[2]) -Base.:(∘)(f::Function,g::Tuple{DistributedCellField,TransientDistributedCellField}) = Operation(f)(g[1],g[2]) - -# Skeleton related -function Base.getproperty(f::TransientDistributedCellField, sym::Symbol) - if sym in (:⁺,:plus,:⁻, :minus) - derivatives = () - cellfield = DistributedCellField(f.cellfield,sym) - for iderivative in f.derivatives - derivatives = (derivatives...,DistributedCellField(iderivative)) - end - return TransientSingleFieldCellField(cellfield,derivatives) - else - return getfield(f, sym) - end -end - -Base.propertynames(x::TransientDistributedCellField, private::Bool=false) = propertynames(x.cellfield, private) - -for op in (:outer,:*,:dot) - @eval begin - ($op)(a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = ($op)(a.cellfield,b) - ($op)(a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = ($op)(a,b.cellfield) - end -end - -Arrays.evaluate!(cache,k::Operation,a::TransientDistributedCellField,b::SkeletonPair{<:DistributedCellField}) = evaluate!(cache,k,a.cellfield,b) - -Arrays.evaluate!(cache,k::Operation,a::SkeletonPair{<:DistributedCellField},b::TransientDistributedCellField) = evaluate!(cache,k,a,b.cellfield) - -CellData.jump(a::TransientDistributedCellField) = jump(a.cellfield) -CellData.mean(a::TransientDistributedCellField) = mean(a.cellfield) diff --git a/src/TransientFESpaces.jl b/src/TransientFESpaces.jl index 6b57a014..980db485 100644 --- a/src/TransientFESpaces.jl +++ b/src/TransientFESpaces.jl @@ -12,47 +12,3 @@ ODEs.∂tt(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.field_f function TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace}) MultiFieldFESpace(spaces) end - -# Functions for transient FE Functions - -function ODEs.allocate_jacobian( - op::ODEs.TransientFEOpFromWeakForm, - t0::Real, - duh::Union{DistributedCellField,DistributedMultiFieldFEFunction}, - cache) - _matdata_jacobians = ODEs.fill_initial_jacobians(op,t0,duh) - matdata = _vcat_distributed_matdata(_matdata_jacobians) - allocate_matrix(op.assem_t,matdata) -end - -function ODEs.jacobians!( - A::AbstractMatrix, - op::ODEs.TransientFEOpFromWeakForm, - t::Real, - xh::TransientDistributedCellField, - γ::Tuple{Vararg{Real}}, - cache) - _matdata_jacobians = ODEs.fill_jacobians(op,t,xh,γ) - matdata = _vcat_distributed_matdata(_matdata_jacobians) - assemble_matrix_add!(A,op.assem_t, matdata) - A -end - -function _vcat_distributed_matdata(_matdata) - term_to_cellmat = map(a->a[1],local_views(_matdata[1])) - term_to_cellidsrows = map(a->a[2],local_views(_matdata[1])) - term_to_cellidscols = map(a->a[3],local_views(_matdata[1])) - for j in 2:length(_matdata) - term_to_cellmat_j = map(a->a[1],local_views(_matdata[j])) - term_to_cellidsrows_j = map(a->a[2],local_views(_matdata[j])) - term_to_cellidscols_j = map(a->a[3],local_views(_matdata[j])) - term_to_cellmat = map((a,b)->vcat(a,b),local_views(term_to_cellmat),local_views(term_to_cellmat_j)) - term_to_cellidsrows = map((a,b)->vcat(a,b),local_views(term_to_cellidsrows),local_views(term_to_cellidsrows_j)) - term_to_cellidscols = map((a,b)->vcat(a,b),local_views(term_to_cellidscols),local_views(term_to_cellidscols_j)) - end - map( (a,b,c) -> (a,b,c), - local_views(term_to_cellmat), - local_views(term_to_cellidsrows), - local_views(term_to_cellidscols) - ) -end From 077d0e3f92eb930d931c81e46eae89fba450952d Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 20 Dec 2023 01:11:17 +1100 Subject: [PATCH 06/25] Completing the transient API --- src/MultiField.jl | 8 +- src/ODEs.jl | 121 ++++++++++++++++-- src/TransientFESpaces.jl | 14 -- ...TransientMultiFieldDistributedCellField.jl | 63 --------- 4 files changed, 113 insertions(+), 93 deletions(-) delete mode 100644 src/TransientFESpaces.jl delete mode 100644 src/TransientMultiFieldDistributedCellField.jl diff --git a/src/MultiField.jl b/src/MultiField.jl index 40f61d51..52b2e0e7 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -107,8 +107,7 @@ function FESpaces.zero_dirichlet_values(f::DistributedMultiFieldFESpace) map(zero_dirichlet_values,f.field_fe_space) end -function FESpaces.FEFunction( - f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false) +function FESpaces.FEFunction(f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false) free_values = change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true) part_fe_fun = map(FEFunction,f.part_fe_space,partition(free_values)) field_fe_fun = DistributedSingleFieldFEFunction[] @@ -121,8 +120,7 @@ function FESpaces.FEFunction( DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) end -function FESpaces.EvaluationFunction( - f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false) +function FESpaces.EvaluationFunction(f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false) free_values = change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true) part_fe_fun = map(EvaluationFunction,f.part_fe_space,partition(free_values)) field_fe_fun = DistributedSingleFieldFEFunction[] @@ -228,7 +226,7 @@ end # DistributedMultiFieldFEBasis -const DistributedMultiFieldFEBasis{A,B<:AbstractArray{<:FEBasis}} = DistributedMultiFieldCellField{A,B} +const DistributedMultiFieldFEBasis{A} = DistributedMultiFieldCellField{A,<:AbstractArray{<:FEBasis}} function FESpaces.get_fe_basis(f::DistributedMultiFieldFESpace) part_mbasis = map(get_fe_basis,f.part_fe_space) diff --git a/src/ODEs.jl b/src/ODEs.jl index d94ddd80..c7e4f4e2 100644 --- a/src/ODEs.jl +++ b/src/ODEs.jl @@ -1,5 +1,62 @@ -# SingleField +# Distributed FESpace + +function Arrays.evaluate!(transient_space::DistributedFESpace, space::DistributedFESpace, t::Real) + map(local_values(transient_space),local_views(space)) do transient_space, space + Arrays.evaluate!(transient_space,space,t) + end + return space +end + +# SingleField FESpace + +const DistributedTransientTrialFESpace = DistributedSingleFieldFESpace{AbstractArray{<:ODEs.AbstractTransientTrialFESpace}} + +function ODEs.TransientTrialFESpace(space::DistributedSingleFieldFESpace,args...) + spaces = map(local_views(space)) do space + ODEs.TransientTrialFESpace(space,args...) + end + gids = get_free_dof_ids(space) + trian = get_triangulation(space) + vector_type = get_vector_type(space) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) +end + +function ODEs.allocate_space(space::DistributedTransientTrialFESpace) + spaces = map(local_views(space)) do space + ODEs.allocate_space(space) + end + gids = get_free_dof_ids(space) + trian = get_triangulation(space) + vector_type = get_vector_type(space) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) +end + +function ODEs.time_derivative(space::DistributedTransientTrialFESpace) + spaces = map(ODEs.time_derivative,local_views(space)) + gids = get_free_dof_ids(space) + trian = get_triangulation(space) + vector_type = get_vector_type(space) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) +end + +for T in [:Real,:Nothing] + @eval begin + function Arrays.evaluate(space::DistributedTransientTrialFESpace, t::$T) + spaces = map(local_views(space)) do space + Arrays.evaluate(space,t) + end + gids = get_free_dof_ids(space) + trian = get_triangulation(space) + vector_type = get_vector_type(space) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) + end + end +end + +# SingleField CellField + +const DistributedTransientSingleFieldCellField{A} = DistributedCellField{A,<:AbstractArray{<:ODEs.TransientSingleFieldCellField}} function ODEs.TransientCellField(f::DistributedCellField,derivatives::Tuple) fields = map(local_views(f),map(local_views,derivatives)...) do f, derivatives... @@ -8,30 +65,72 @@ function ODEs.TransientCellField(f::DistributedCellField,derivatives::Tuple) DistributedCellField(fields,get_triangulation(f),f.metadata) end -function ODEs.time_derivative(f::DistributedCellField) +function ODEs.time_derivative(f::DistributedTransientSingleFieldCellField) fields = map(local_views(f)) do field ODEs.time_derivative(field) end DistributedCellField(fields,get_triangulation(f)) end -# MultiField +# MultiField FESpace + +const DistributedTransientMultiFieldFESpace{MS,A} = + DistributedMultiFieldFESpace{MS,A,<:AbstractArray{<:ODEs.TransientMultiFieldFESpace}} + +function ODEs.TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace}) + MultiFieldFESpace(spaces) +end + +function ODEs.allocate_space(space::DistributedTransientMultiFieldFESpace{MS}) where MS + field_fe_spaces = map(ODEs.allocate_space,space.field_fe_spaces) + spaces = to_parray_of_arrays(map(local_views,field_fe_spaces)) + part_fe_spaces = map(s -> MultiFieldFESpace(s;style=MS()),spaces) + gids = get_free_dof_ids(space) + vector_type = get_vector_type(space) + DistributedMultiFieldFESpace(field_fe_spaces,part_fe_spaces,gids,vector_type) +end + +function ODEs.time_derivative(f::DistributedTransientMultiFieldFESpace{MS}) where MS + field_fe_spaces = map(ODEs.time_derivative,f.field_fe_spaces) + spaces = to_parray_of_arrays(map(local_views,field_fe_spaces)) + part_fe_spaces = map(s -> MultiFieldFESpace(s;style=MS()),spaces) + gids = get_free_dof_ids(f) + vector_type = get_vector_type(f) + DistributedMultiFieldFESpace(field_fe_spaces,part_fe_spaces,gids,vector_type) +end + +for T in [:Real,:Nothing] + @eval begin + function Arrays.evaluate(space::DistributedTransientMultiFieldFESpace{MS}, t::$T) where MS + field_fe_spaces = map(s->Arrays.evaluate(s,t),space.field_fe_spaces) + spaces = to_parray_of_arrays(map(local_views,field_fe_spaces)) + part_fe_spaces = map(s -> MultiFieldFESpace(s;style=MS()),spaces) + gids = get_free_dof_ids(space) + vector_type = get_vector_type(space) + DistributedMultiFieldFESpace(field_fe_spaces,part_fe_spaces,gids,vector_type) + end + end +end -function ODEs.TransientCellField(f::DistributedMultiFieldFEFunction,derivatives::Tuple) +# MultiField CellField + +const DistributedTransientMultiFieldCellField{A} = + DistributedMultiFieldCellField{A,<:AbstractArray{<:ODEs.TransientMultiFieldCellField}} + +function ODEs.TransientCellField(f::DistributedMultiFieldCellField,derivatives::Tuple) field_fe_fun = map(1:num_fields(f)) do i f_i = f[i] df_i = Tuple(map(df -> df[i],derivatives)) ODEs.TransientCellField(f_i,df_i) end - part_fe_fun = map(local_views(f),map(local_views,derivatives)...) do f, derivatives... - ODEs.TransientCellField(f,Tuple(derivatives)) - end - free_values = get_free_dof_values(f) - DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) + fields = to_parray_of_arrays(map(local_views,field_fe_fun)) + part_fe_fun = map(MultiField.MultiFieldCellField,fields) + DistributedMultiFieldCellField(field_fe_fun,part_fe_fun,f.metadata) end -function ODEs.time_derivative(f::DistributedMultiFieldCellField) +function ODEs.time_derivative(f::DistributedTransientMultiFieldCellField) field_fe_fun = map(ODEs.time_derivative,f.field_fe_fun) - part_fe_fun = map(ODEs.time_derivative,f.part_fe_fun) + fields = to_parray_of_arrays(map(local_views,field_fe_fun)) + part_fe_fun = map(MultiField.MultiFieldCellField,fields) DistributedMultiFieldCellField(field_fe_fun,part_fe_fun) end diff --git a/src/TransientFESpaces.jl b/src/TransientFESpaces.jl deleted file mode 100644 index 980db485..00000000 --- a/src/TransientFESpaces.jl +++ /dev/null @@ -1,14 +0,0 @@ -# Functions for transient FE spaces - -Fields.evaluate(U::DistributedSingleFieldFESpace,t::Nothing) = U - -(U::DistributedSingleFieldFESpace)(t) = U - -ODEs.∂t(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U) -ODEs.∂t(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂t.(U.field_fe_space)) -ODEs.∂tt(U::DistributedSingleFieldFESpace) = HomogeneousTrialFESpace(U) -ODEs.∂tt(U::DistributedMultiFieldFESpace) = MultiFieldFESpace(∂tt.(U.field_fe_spaces)) - -function TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace}) - MultiFieldFESpace(spaces) -end diff --git a/src/TransientMultiFieldDistributedCellField.jl b/src/TransientMultiFieldDistributedCellField.jl deleted file mode 100644 index 3f4fd47f..00000000 --- a/src/TransientMultiFieldDistributedCellField.jl +++ /dev/null @@ -1,63 +0,0 @@ -# Transient MultiField -struct TransientMultiFieldDistributedCellField{A} <: TransientDistributedCellField - cellfield::A - derivatives::Tuple - transient_single_fields::Vector{<:TransientDistributedCellField} # used to iterate -end - -local_views(f::TransientMultiFieldDistributedCellField) = local_views(f.cellfield) - -# Constructors -function ODEs.TransientCellField(multi_field::DistributedMultiFieldFEFunction,derivatives::Tuple) - transient_single_fields = _to_transient_single_distributed_fields(multi_field,derivatives) - TransientMultiFieldDistributedCellField(multi_field,derivatives,transient_single_fields) -end - -# Get single index -function Base.getindex(f::TransientMultiFieldDistributedCellField,ifield::Integer) - single_field = f.cellfield[ifield] - single_derivatives = () - for ifield_derivatives in f.derivatives - single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield)) - end - TransientSingleFieldDistributedCellField(single_field,single_derivatives) -end - -# Get multiple indices -function Base.getindex(f::TransientMultiFieldDistributedCellField,indices::Vector{<:Int}) - cellfield = DistributedMultiFieldCellField(f.cellfield[indices],DomainStyle(f.cellfield)) - derivatives = () - for derivative in f.derivatives - derivatives = (derivatives...,DistributedMultiFieldCellField(derivative[indices],DomainStyle(derivative))) - end - transient_single_fields = _to_transient_single_distributed_fields(cellfield,derivatives) - TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_fields) -end - -function _to_transient_single_distributed_fields(multi_field,derivatives) - transient_single_fields = TransientDistributedCellField[] - for ifield in 1:num_fields(multi_field) - single_field = multi_field[ifield] - single_derivatives = () - for ifield_derivatives in derivatives - single_derivatives = (single_derivatives...,getindex(ifield_derivatives,ifield)) - end - transient_single_field = TransientSingleFieldDistributedCellField(single_field,single_derivatives) - push!(transient_single_fields,transient_single_field) - end - transient_single_fields -end - -# Iterate functions -Base.iterate(f::TransientMultiFieldDistributedCellField) = iterate(f.transient_single_fields) -Base.iterate(f::TransientMultiFieldDistributedCellField,state) = iterate(f.transient_single_fields,state) - -# Time derivative -function ODEs.∂t(f::TransientMultiFieldDistributedCellField) - cellfield, derivatives = first_and_tail(f.derivatives) - transient_single_field_derivatives = TransientDistributedCellField[] - for transient_single_field in f.transient_single_fields - push!(transient_single_field_derivatives,∂t(transient_single_field)) - end - TransientMultiFieldDistributedCellField(cellfield,derivatives,transient_single_field_derivatives) -end From 514ae587b5c4b4428bc3c0b5c1e02c258be57024 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 20 Dec 2023 01:19:00 +1100 Subject: [PATCH 07/25] Minor --- src/ODEs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ODEs.jl b/src/ODEs.jl index c7e4f4e2..6552177b 100644 --- a/src/ODEs.jl +++ b/src/ODEs.jl @@ -1,5 +1,5 @@ -# Distributed FESpace +# Distributed FESpace commons function Arrays.evaluate!(transient_space::DistributedFESpace, space::DistributedFESpace, t::Real) map(local_values(transient_space),local_views(space)) do transient_space, space From 7edb9e5bd3f1e281dbb54c7c27b2eafb2b9df687 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 20 Dec 2023 09:33:28 +1100 Subject: [PATCH 08/25] Fixes for transient multifield --- src/ODEs.jl | 14 +++++++++++--- test/TransientDistributedCellFieldTests.jl | 9 +++------ ...sientMultiFieldDistributedCellFieldTests.jl | 18 ++++++------------ 3 files changed, 20 insertions(+), 21 deletions(-) diff --git a/src/ODEs.jl b/src/ODEs.jl index 6552177b..866ab8e6 100644 --- a/src/ODEs.jl +++ b/src/ODEs.jl @@ -56,7 +56,7 @@ end # SingleField CellField -const DistributedTransientSingleFieldCellField{A} = DistributedCellField{A,<:AbstractArray{<:ODEs.TransientSingleFieldCellField}} +const DistributedTransientSingleFieldCellField = DistributedCellField{<:AbstractArray{<:ODEs.TransientSingleFieldCellField}} function ODEs.TransientCellField(f::DistributedCellField,derivatives::Tuple) fields = map(local_views(f),map(local_views,derivatives)...) do f, derivatives... @@ -124,13 +124,21 @@ function ODEs.TransientCellField(f::DistributedMultiFieldCellField,derivatives:: ODEs.TransientCellField(f_i,df_i) end fields = to_parray_of_arrays(map(local_views,field_fe_fun)) - part_fe_fun = map(MultiField.MultiFieldCellField,fields) + part_fe_fun = map(ODEs.TransientMultiFieldCellField,fields) DistributedMultiFieldCellField(field_fe_fun,part_fe_fun,f.metadata) end +function ODEs.TransientMultiFieldCellField(fields::AbstractVector{<:ODEs.TransientSingleFieldCellField}) + cellfield = MultiFieldCellField(map(f -> f.cellfield,fields)) + n_derivatives = length(first(fields).derivatives) + @check all(map(f -> length(f.derivatives) == n_derivatives,fields)) + derivatives = Tuple(map(i -> MultiFieldCellField(map(f -> f.derivatives[i],fields)),1:n_derivatives)) + TransientMultiFieldCellField(cellfield,derivatives,fields) +end + function ODEs.time_derivative(f::DistributedTransientMultiFieldCellField) field_fe_fun = map(ODEs.time_derivative,f.field_fe_fun) fields = to_parray_of_arrays(map(local_views,field_fe_fun)) - part_fe_fun = map(MultiField.MultiFieldCellField,fields) + part_fe_fun = map(ODEs.TransientMultiFieldCellField,fields) DistributedMultiFieldCellField(field_fe_fun,part_fe_fun) end diff --git a/test/TransientDistributedCellFieldTests.jl b/test/TransientDistributedCellFieldTests.jl index aed283d1..b6328677 100644 --- a/test/TransientDistributedCellFieldTests.jl +++ b/test/TransientDistributedCellFieldTests.jl @@ -27,16 +27,13 @@ function main(distribute,parts) @test isa(dda(0),GridapDistributed.DistributedCellField) b(t) = TransientCellField(a(t),(da(t),dda(t))) - @test isa(b(0),GridapDistributed.TransientDistributedCellField) - @test isa(b(0),GridapDistributed.TransientSingleFieldDistributedCellField) + @test isa(b(0),GridapDistributed.DistributedTransientSingleFieldCellField) db(t) = ∂t(b(t)) - @test isa(db(0),GridapDistributed.TransientDistributedCellField) - @test isa(db(0),GridapDistributed.TransientSingleFieldDistributedCellField) + @test isa(db(0),GridapDistributed.DistributedTransientSingleFieldCellField) ddb(t) = ∂t(db(t)) - @test isa(ddb(0),GridapDistributed.TransientDistributedCellField) - @test isa(ddb(0),GridapDistributed.TransientSingleFieldDistributedCellField) + @test isa(ddb(0),GridapDistributed.DistributedTransientSingleFieldCellField) @test (∑(∫(a(0.5))dΩ)) ≈ 0.25 @test (∑(∫(da(0.5))dΩ)) ≈ 1.0 diff --git a/test/TransientMultiFieldDistributedCellFieldTests.jl b/test/TransientMultiFieldDistributedCellFieldTests.jl index 76b66a76..ca83333e 100644 --- a/test/TransientMultiFieldDistributedCellFieldTests.jl +++ b/test/TransientMultiFieldDistributedCellFieldTests.jl @@ -36,28 +36,22 @@ function main(distribute,parts) @test isa(dda(0),GridapDistributed.DistributedMultiFieldFEFunction) b(t) = TransientCellField(a(t),(da(t),dda(t))) - @test isa(b(0),GridapDistributed.TransientDistributedCellField) - @test isa(b(0),GridapDistributed.TransientMultiFieldDistributedCellField) + @test isa(b(0),GridapDistributed.DistributedTransientMultiFieldCellField) db(t) = ∂t(b(t)) - @test isa(db(0),GridapDistributed.TransientDistributedCellField) - @test isa(db(0),GridapDistributed.TransientMultiFieldDistributedCellField) + @test isa(db(0),GridapDistributed.DistributedTransientMultiFieldCellField) ddb(t) = ∂t(db(t)) - @test isa(ddb(0),GridapDistributed.TransientDistributedCellField) - @test isa(ddb(0),GridapDistributed.TransientMultiFieldDistributedCellField) + @test isa(ddb(0),GridapDistributed.DistributedTransientMultiFieldCellField) b1(t) = b(t)[1] - @test isa(b1(0),GridapDistributed.TransientDistributedCellField) - @test isa(b1(0),GridapDistributed.TransientSingleFieldDistributedCellField) + @test isa(b1(0),GridapDistributed.DistributedTransientSingleFieldCellField) db1(t) = ∂t(b1(t)) - @test isa(db1(0),GridapDistributed.TransientDistributedCellField) - @test isa(db1(0),GridapDistributed.TransientSingleFieldDistributedCellField) + @test isa(db1(0),GridapDistributed.DistributedTransientSingleFieldCellField) ddb1(t) = ∂t(db1(t)) - @test isa(ddb1(0),GridapDistributed.TransientDistributedCellField) - @test isa(ddb1(0),GridapDistributed.TransientSingleFieldDistributedCellField) + @test isa(ddb1(0),GridapDistributed.DistributedTransientSingleFieldCellField) @test (∑(∫(b(0.5)[1])dΩ)) == (∑(∫(b1(0.5))dΩ)) @test (∑(∫(db(0.5)[1])dΩ)) == (∑(∫(db1(0.5))dΩ)) From 85e9765d36ce2ca88c8f766eb0bacca4bf7d0116 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 20 Dec 2023 16:59:56 +1100 Subject: [PATCH 09/25] First working version --- src/Algebra.jl | 40 ++++++++++++++++++ src/CellData.jl | 11 +++++ src/MultiField.jl | 1 + src/ODEs.jl | 75 +++++++++++++++++++++------------ test/HeatEquationTests.jl | 28 ++++++------ test/StokesOpenBoundaryTests.jl | 31 +++++--------- 6 files changed, 125 insertions(+), 61 deletions(-) diff --git a/src/Algebra.jl b/src/Algebra.jl index 45244eca..33d7c293 100644 --- a/src/Algebra.jl +++ b/src/Algebra.jl @@ -29,6 +29,46 @@ function Algebra.allocate_in_domain(matrix::BlockPMatrix) allocate_in_domain(BlockPVector{V},matrix) end +# PartitionedArrays extras + +function LinearAlgebra.axpy!(α,x::PVector,y::PVector) + @check partition(axes(x,1)) === partition(axes(y,1)) + map(partition(x),partition(y)) do x,y + LinearAlgebra.axpy!(α,x,y) + end + return y +end + +function LinearAlgebra.axpy!(α,x::BlockPVector,y::BlockPVector) + map(blocks(x),blocks(y)) do x,y + LinearAlgebra.axpy!(α,x,y) + end + return y +end + +function Algebra.axpy_entries!( + α::Number, A::PSparseMatrix, B::PSparseMatrix; + check::Bool=true +) +# We should definitely check here that the index partitions are the same. +# However: Because the different matrices are assembled separately, the objects are not the +# same (i.e can't use ===). Checking the index partitions would then be costly... + map(partition(A),partition(B)) do A, B + Algebra.axpy_entries!(α,A,B;check) + end + return B +end + +function Algebra.axpy_entries!( + α::Number, A::BlockPMatrix, B::BlockPMatrix; + check::Bool=true +) + map(blocks(A),blocks(B)) do A, B + Algebra.axpy_entries!(α,A,B;check) + end + return B +end + # This might go to Gridap in the future. We keep it here for the moment. function change_axes(a::Algebra.ArrayCounter,axes) @notimplemented diff --git a/src/CellData.jl b/src/CellData.jl index 68983440..0814b398 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -221,6 +221,8 @@ struct DistributedDomainContribution{A<:AbstractArray{<:DomainContribution}} <: contribs::A end +CellData.num_domains(a::DistributedDomainContribution) = CellData.num_domains(getany(local_views(a))) + local_views(a::DistributedDomainContribution) = a.contribs function Base.getindex(c::DistributedDomainContribution,t::DistributedTriangulation) @@ -277,6 +279,15 @@ end (*)(a::DistributedDomainContribution,b::Number) = b*a +# Jordi: This is ugly, but it is useful to re-use code from Gridap: +# A lot of the time, we create an empty DomainContribution and then add to it. +# By dispatching here, this kind of code works verbatim for GridapDistributed. +# We could eventually replace this with an EmptyDomainContribution type. +function (+)(a::CellData.DomainContribution,b::DistributedDomainContribution) + @assert iszero(CellData.num_domains(a)) + return b +end + # Triangulation related function CellData.get_cell_points(a::DistributedTriangulation) diff --git a/src/MultiField.jl b/src/MultiField.jl index 52b2e0e7..a2a4ae1e 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -78,6 +78,7 @@ function CellData.get_triangulation(a::DistributedMultiFieldFESpace) end MultiField.MultiFieldStyle(::Type{<:DistributedMultiFieldFESpace{MS}}) where MS = MS() +MultiField.MultiFieldStyle(a::DistributedMultiFieldFESpace) = MultiField.MultiFieldStyle(typeof(a)) local_views(a::DistributedMultiFieldFESpace) = a.part_fe_space MultiField.num_fields(m::DistributedMultiFieldFESpace) = length(m.field_fe_space) diff --git a/src/ODEs.jl b/src/ODEs.jl index 866ab8e6..1073ea4b 100644 --- a/src/ODEs.jl +++ b/src/ODEs.jl @@ -2,19 +2,29 @@ # Distributed FESpace commons function Arrays.evaluate!(transient_space::DistributedFESpace, space::DistributedFESpace, t::Real) - map(local_values(transient_space),local_views(space)) do transient_space, space + map(local_views(transient_space),local_views(space)) do transient_space, space Arrays.evaluate!(transient_space,space,t) end - return space + return transient_space end # SingleField FESpace -const DistributedTransientTrialFESpace = DistributedSingleFieldFESpace{AbstractArray{<:ODEs.AbstractTransientTrialFESpace}} +const DistributedTransientTrialFESpace = DistributedSingleFieldFESpace{<:AbstractArray{<:ODEs.TransientTrialFESpace}} -function ODEs.TransientTrialFESpace(space::DistributedSingleFieldFESpace,args...) +function ODEs.TransientTrialFESpace(space::DistributedSingleFieldFESpace,transient_dirichlet::Union{Function, AbstractVector{<:Function}}) spaces = map(local_views(space)) do space - ODEs.TransientTrialFESpace(space,args...) + ODEs.TransientTrialFESpace(space,transient_dirichlet) + end + gids = get_free_dof_ids(space) + trian = get_triangulation(space) + vector_type = get_vector_type(space) + DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) +end + +function ODEs.TransientTrialFESpace(space::DistributedSingleFieldFESpace) + spaces = map(local_views(space)) do space + ODEs.TransientTrialFESpace(space) end gids = get_free_dof_ids(space) trian = get_triangulation(space) @@ -32,7 +42,7 @@ function ODEs.allocate_space(space::DistributedTransientTrialFESpace) DistributedSingleFieldFESpace(spaces,gids,trian,vector_type) end -function ODEs.time_derivative(space::DistributedTransientTrialFESpace) +function ODEs.time_derivative(space::DistributedSingleFieldFESpace) spaces = map(ODEs.time_derivative,local_views(space)) gids = get_free_dof_ids(space) trian = get_triangulation(space) @@ -74,40 +84,49 @@ end # MultiField FESpace -const DistributedTransientMultiFieldFESpace{MS,A} = - DistributedMultiFieldFESpace{MS,A,<:AbstractArray{<:ODEs.TransientMultiFieldFESpace}} - -function ODEs.TransientMultiFieldFESpace(spaces::Vector{<:DistributedSingleFieldFESpace}) - MultiFieldFESpace(spaces) +function ODEs.has_transient(space::DistributedMultiFieldFESpace) + getany(map(ODEs.has_transient,local_views(space))) end -function ODEs.allocate_space(space::DistributedTransientMultiFieldFESpace{MS}) where MS - field_fe_spaces = map(ODEs.allocate_space,space.field_fe_spaces) - spaces = to_parray_of_arrays(map(local_views,field_fe_spaces)) - part_fe_spaces = map(s -> MultiFieldFESpace(s;style=MS()),spaces) +function ODEs.allocate_space(space::DistributedMultiFieldFESpace) + if !ODEs.has_transient(space) + return space + end + field_fe_space = map(ODEs.allocate_space,space.field_fe_space) + style = MultiFieldStyle(space) + spaces = to_parray_of_arrays(map(local_views,field_fe_space)) + part_fe_spaces = map(s -> MultiFieldFESpace(s;style),spaces) gids = get_free_dof_ids(space) vector_type = get_vector_type(space) - DistributedMultiFieldFESpace(field_fe_spaces,part_fe_spaces,gids,vector_type) + DistributedMultiFieldFESpace(field_fe_space,part_fe_spaces,gids,vector_type) end -function ODEs.time_derivative(f::DistributedTransientMultiFieldFESpace{MS}) where MS - field_fe_spaces = map(ODEs.time_derivative,f.field_fe_spaces) - spaces = to_parray_of_arrays(map(local_views,field_fe_spaces)) - part_fe_spaces = map(s -> MultiFieldFESpace(s;style=MS()),spaces) - gids = get_free_dof_ids(f) - vector_type = get_vector_type(f) - DistributedMultiFieldFESpace(field_fe_spaces,part_fe_spaces,gids,vector_type) +function ODEs.time_derivative(space::DistributedMultiFieldFESpace) + if !ODEs.has_transient(space) + return space + end + field_fe_space = map(ODEs.time_derivative,space.field_fe_space) + style = MultiFieldStyle(space) + spaces = to_parray_of_arrays(map(local_views,field_fe_space)) + part_fe_spaces = map(s -> MultiFieldFESpace(s;style),spaces) + gids = get_free_dof_ids(space) + vector_type = get_vector_type(space) + DistributedMultiFieldFESpace(field_fe_space,part_fe_spaces,gids,vector_type) end for T in [:Real,:Nothing] @eval begin - function Arrays.evaluate(space::DistributedTransientMultiFieldFESpace{MS}, t::$T) where MS - field_fe_spaces = map(s->Arrays.evaluate(s,t),space.field_fe_spaces) - spaces = to_parray_of_arrays(map(local_views,field_fe_spaces)) - part_fe_spaces = map(s -> MultiFieldFESpace(s;style=MS()),spaces) + function Arrays.evaluate(space::DistributedMultiFieldFESpace, t::$T) + if !ODEs.has_transient(space) + return space + end + field_fe_space = map(s->Arrays.evaluate(s,t),space.field_fe_space) + style = MultiFieldStyle(space) + spaces = to_parray_of_arrays(map(local_views,field_fe_space)) + part_fe_spaces = map(s -> MultiFieldFESpace(s;style),spaces) gids = get_free_dof_ids(space) vector_type = get_vector_type(space) - DistributedMultiFieldFESpace(field_fe_spaces,part_fe_spaces,gids,vector_type) + DistributedMultiFieldFESpace(field_fe_space,part_fe_spaces,gids,vector_type) end end end diff --git a/test/HeatEquationTests.jl b/test/HeatEquationTests.jl index 403d8ead..80fc9fec 100644 --- a/test/HeatEquationTests.jl +++ b/test/HeatEquationTests.jl @@ -1,7 +1,7 @@ module HeatEquationTests using Gridap -using Gridap.ODEs +using Gridap.ODEs, Gridap.Algebra using GridapDistributed using PartitionedArrays using Test @@ -34,41 +34,43 @@ function main(distribute,parts) dΩ = Measure(Ω,degree) # - m(u,v) = ∫(u*v)dΩ - a(u,v) = ∫(∇(v)⋅∇(u))dΩ + m(t,u,v) = ∫(u*v)dΩ + a(t,u,v) = ∫(∇(v)⋅∇(u))dΩ b(t,v) = ∫(v*f(t))dΩ - res(t,u,v) = a(u,v) + m(∂t(u),v) - b(t,v) - jac(t,u,du,v) = a(du,v) - jac_t(t,u,dut,v) = m(dut,v) + res(t,u,v) = a(t,u,v) + m(t,∂t(u),v) - b(t,v) + jac(t,u,du,v) = a(t,du,v) + jac_t(t,u,dut,v) = m(t,dut,v) op = TransientFEOperator(res,jac,jac_t,U,V0) - op_constant = TransientConstantMatrixFEOperator(m,a,b,U,V0) + op_constant = TransientLinearFEOperator(a,m,(t,v) -> (-1)*b(t,v),U,V0,constant_forms=(true,true)) t0 = 0.0 tF = 1.0 dt = 0.1 U0 = U(0.0) - uh0 = interpolate_everywhere(u(0.0),U0) + uh0 = interpolate_everywhere(u(0.0),U0); ls = LUSolver() - ode_solver = ThetaMethod(ls,dt,θ) + linear_ode_solver = ThetaMethod(ls,dt,θ) + sol_t_const = solve(linear_ode_solver,op_constant,t0,tF,uh0) - sol_t = solve(ode_solver,op,uh0,t0,tF) - sol_t_const = solve(ode_solver,op_constant,uh0,t0,tF) + nls = NewtonRaphsonSolver(ls,1.0e-6,10) + nonlinear_ode_solver = ThetaMethod(nls,dt,θ) + sol_t = solve(nonlinear_ode_solver,op,t0,tF,uh0) l2(w) = w*w tol = 1.0e-6 - for (uh_tn, tn) in sol_t + for (tn, uh_tn) in sol_t e = u(tn) - uh_tn el2 = sqrt(sum( ∫(l2(e))dΩ )) @test el2 < tol end - for (uh_tn, tn) in sol_t_const + for (tn, uh_tn) in sol_t_const e = u(tn) - uh_tn el2 = sqrt(sum( ∫(l2(e))dΩ )) @test el2 < tol diff --git a/test/StokesOpenBoundaryTests.jl b/test/StokesOpenBoundaryTests.jl index 93bf733a..d3a61b30 100644 --- a/test/StokesOpenBoundaryTests.jl +++ b/test/StokesOpenBoundaryTests.jl @@ -59,20 +59,16 @@ function main(distribute,parts) dΓ = Measure(Γ,degree) # - a(u,v) = ∫(∇(u)⊙∇(v))dΩ - b((v,q),t) = ∫(v⋅f(t))dΩ + ∫(q*g(t))dΩ + ∫(v⋅h(t))dΓ - m(ut,v) = ∫(ut⋅v)dΩ + a(t,u,v) = ∫(∇(u)⊙∇(v))dΩ + b(t,(v,q)) = ∫(v⋅f(t))dΩ + ∫(q*g(t))dΩ + ∫(v⋅h(t))dΓ + m(t,ut,v) = ∫(ut⋅v)dΩ X = TransientMultiFieldFESpace([U,P]) Y = MultiFieldFESpace([V0,Q]) - res(t,(u,p),(v,q)) = a(u,v) + m(∂t(u),v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b((v,q),t) - jac(t,(u,p),(du,dp),(v,q)) = a(du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ - jac_t(t,(u,p),(dut,dpt),(v,q)) = m(dut,v) - - b((v,q)) = b((v,q),0.0) - - mat((du1,du2),(v1,v2)) = a(du1,v1)+a(du2,v2) + res(t,(u,p),(v,q)) = a(t,u,v) + m(t,∂t(u),v) - ∫((∇⋅v)*p)dΩ + ∫(q*(∇⋅u))dΩ - b(t,(v,q)) + jac(t,(u,p),(du,dp),(v,q)) = a(t,du,v) - ∫((∇⋅v)*dp)dΩ + ∫(q*(∇⋅du))dΩ + jac_t(t,(u,p),(dut,dpt),(v,q)) = m(t,dut,v) U0 = U(0.0) P0 = P(0.0) @@ -87,23 +83,18 @@ function main(distribute,parts) tF = 1.0 dt = 0.1 - ls = LUSolver() - ode_solver = ThetaMethod(ls,dt,θ) + ls = LUSolver() + nls = NewtonRaphsonSolver(ls,1.0e-6,10) + ode_solver = ThetaMethod(nls,dt,θ) sol_t = solve(ode_solver,op,t0,tF,xh0) l2(w) = w⋅w - - tol = 1.0e-6 - _t_n = t0 - - result = Base.iterate(sol_t) - - for (xh_tn, tn) in sol_t + for (tn, xh_tn) in sol_t uh_tn = xh_tn[1] ph_tn = xh_tn[2] - writevtk(Ω,"output/tmp_stokes_OB_sol_$tn.vtu",cellfields=["u"=>uh_tn,"p"=>ph_tn]) + #writevtk(Ω,"output/tmp_stokes_OB_sol_$tn.vtu",cellfields=["u"=>uh_tn,"p"=>ph_tn]) e = u(tn) - uh_tn el2 = sqrt(sum( ∫(l2(e))dΩ )) e = p(tn) - ph_tn From 0372ca68a0e33efb71713421b8f3d266e2104647 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 20 Mar 2024 09:37:23 +1100 Subject: [PATCH 10/25] Added some missing API for MultiField --- src/MultiField.jl | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/src/MultiField.jl b/src/MultiField.jl index a2a4ae1e..56114b13 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -121,6 +121,23 @@ function FESpaces.FEFunction(f::DistributedMultiFieldFESpace,x::AbstractVector,i DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) end +function FESpaces.FEFunction( + f::DistributedMultiFieldFESpace,x::AbstractVector, + dirichlet_values::AbstractArray{<:AbstractVector},isconsistent=false + ) + free_values = GridapDistributed.change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true) + part_fe_fun = map(FEFunction,f.part_fe_space,partition(free_values)) + field_fe_fun = DistributedSingleFieldFEFunction[] + for i in 1:num_fields(f) + free_values_i = restrict_to_field(f,free_values,i) + dirichlet_values_i = dirichlet_values[i] + fe_space_i = f.field_fe_space[i] + fe_fun_i = FEFunction(fe_space_i,free_values_i,dirichlet_values_i,true) + push!(field_fe_fun,fe_fun_i) + end + DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) +end + function FESpaces.EvaluationFunction(f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false) free_values = change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true) part_fe_fun = map(EvaluationFunction,f.part_fe_space,partition(free_values)) @@ -153,6 +170,24 @@ function FESpaces.interpolate!(objects,free_values::AbstractVector,fe::Distribut DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) end +function Gridap.FESpaces.interpolate!( + objects::Union{<:DistributedMultiFieldCellField,<:DistributedCellField}, + free_values::AbstractVector, + fe::DistributedMultiFieldFESpace +) + part_fe_fun = map(local_views(objects),partition(free_values),local_views(fe)) do objects,x,f + interpolate!(objects,x,f) + end + field_fe_fun = GridapDistributed.DistributedSingleFieldFEFunction[] + for i in 1:num_fields(fe) + free_values_i = Gridap.MultiField.restrict_to_field(fe,free_values,i) + fe_space_i = fe.field_fe_space[i] + fe_fun_i = FEFunction(fe_space_i,free_values_i) + push!(field_fe_fun,fe_fun_i) + end + GridapDistributed.DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) +end + function FESpaces.interpolate_everywhere(objects,fe::DistributedMultiFieldFESpace) free_values = zero_free_values(fe) part_fe_fun = map(partition(free_values),local_views(fe)) do x,f From 40a7fd2651fa097fb381c11128c778edd176ba04 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 21 Mar 2024 14:34:14 +1100 Subject: [PATCH 11/25] Fixed bug in PArtitionedArrays copy function --- src/Algebra.jl | 24 +++++++++++++++++------- test/HeatEquationTests.jl | 20 +++++++------------- 2 files changed, 24 insertions(+), 20 deletions(-) diff --git a/src/Algebra.jl b/src/Algebra.jl index b191e32c..47d6ca40 100644 --- a/src/Algebra.jl +++ b/src/Algebra.jl @@ -29,6 +29,14 @@ function Algebra.allocate_in_domain(matrix::BlockPMatrix) allocate_in_domain(BlockPVector{V},matrix) end +# PSparseMatrix copy + +function Base.copy(a::PSparseMatrix) + mats = map(copy,partition(a)) + cache = map(PartitionedArrays.copy_cache,a.cache) + return PSparseMatrix(mats,partition(axes(a,1)),partition(axes(a,2)),cache) +end + # PartitionedArrays extras function LinearAlgebra.axpy!(α,x::PVector,y::PVector) @@ -36,6 +44,7 @@ function LinearAlgebra.axpy!(α,x::PVector,y::PVector) map(partition(x),partition(y)) do x,y LinearAlgebra.axpy!(α,x,y) end + consistent!(y) |> wait return y end @@ -53,6 +62,8 @@ function Algebra.axpy_entries!( # We should definitely check here that the index partitions are the same. # However: Because the different matrices are assembled separately, the objects are not the # same (i.e can't use ===). Checking the index partitions would then be costly... + @assert reduce(&,map(PartitionedArrays.matching_local_indices,partition(axes(A,1)),partition(axes(B,1)))) + @assert reduce(&,map(PartitionedArrays.matching_local_indices,partition(axes(A,2)),partition(axes(B,2)))) map(partition(A),partition(B)) do A, B Algebra.axpy_entries!(α,A,B;check) end @@ -835,12 +846,10 @@ function local_views(a::PVectorAllocationTrackTouchedAndValues) end function _setup_prange_from_pvector_allocation(a::PVectorAllocationTrackTouchedAndValues) - test_dofs_prange = a.test_dofs_gids_prange # dof ids of the test space - ngrdofs = length(test_dofs_prange) - + # Find the ghost rows allocations = local_views(a.allocations) - indices = partition(test_dofs_prange) + indices = partition(a.test_dofs_gids_prange) I_ghost_lids_to_dofs_ghost_lids = map(allocations, indices) do allocation, indices dofs_lids_touched = findall(allocation.touched) loc_to_gho = local_to_ghost(indices) @@ -857,10 +866,11 @@ function _setup_prange_from_pvector_allocation(a::PVectorAllocationTrackTouchedA I_ghost_lids end - gids_ghost_to_global, gids_ghost_to_owner = map( + ghost_to_global, ghost_to_owner = map( find_gid_and_owner,I_ghost_lids_to_dofs_ghost_lids,indices) |> tuple_of_arrays - _setup_prange_impl_(ngrdofs,indices,gids_ghost_to_global,gids_ghost_to_owner) + ngids = length(a.test_dofs_gids_prange) + _setup_prange_impl_(ngids,indices,ghost_to_global,ghost_to_owner) end function Algebra.create_from_nz(a::PVectorAllocationTrackTouchedAndValues) @@ -888,7 +898,7 @@ function find_gid_and_owner(ighost_to_jghost,jindices) jghost_to_local = ghost_to_local(jindices) jlocal_to_global = local_to_global(jindices) jlocal_to_owner = local_to_owner(jindices) - ighost_to_jlocal = view(jghost_to_local,ighost_to_jghost) + ighost_to_jlocal = sort(view(jghost_to_local,ighost_to_jghost)) ighost_to_global = jlocal_to_global[ighost_to_jlocal] ighost_to_owner = jlocal_to_owner[ighost_to_jlocal] diff --git a/test/HeatEquationTests.jl b/test/HeatEquationTests.jl index 80fc9fec..6bda0178 100644 --- a/test/HeatEquationTests.jl +++ b/test/HeatEquationTests.jl @@ -15,25 +15,17 @@ function main(distribute,parts) f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x) domain = (0,1,0,1) - partition = (4,4) - model = CartesianDiscreteModel(ranks,parts,domain,partition) + model = CartesianDiscreteModel(ranks,parts,domain,(4,4)) order = 2 - reffe = ReferenceFE(lagrangian,Float64,order) - V0 = FESpace( - model, - reffe, - conformity=:H1, - dirichlet_tags="boundary" - ) + V0 = FESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary") U = TransientTrialFESpace(V0,u) Ω = Triangulation(model) degree = 2*order dΩ = Measure(Ω,degree) - # m(t,u,v) = ∫(u*v)dΩ a(t,u,v) = ∫(∇(v)⋅∇(u))dΩ b(t,v) = ∫(v*f(t))dΩ @@ -43,7 +35,11 @@ function main(distribute,parts) jac_t(t,u,dut,v) = m(t,dut,v) op = TransientFEOperator(res,jac,jac_t,U,V0) - op_constant = TransientLinearFEOperator(a,m,(t,v) -> (-1)*b(t,v),U,V0,constant_forms=(true,true)) + + assembler = SparseMatrixAssembler(U,V0,SubAssembledRows()) + op_constant = TransientLinearFEOperator( + (a,m),(t,v) -> (-1)*b(t,v),U,V0,constant_forms=(true,true),assembler=assembler + ) t0 = 0.0 tF = 1.0 @@ -61,9 +57,7 @@ function main(distribute,parts) sol_t = solve(nonlinear_ode_solver,op,t0,tF,uh0) l2(w) = w*w - tol = 1.0e-6 - for (tn, uh_tn) in sol_t e = u(tn) - uh_tn el2 = sqrt(sum( ∫(l2(e))dΩ )) From ae298c86d1fcf84a1797278c5e4c532fe7b1433a Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 21 Mar 2024 17:44:35 +1100 Subject: [PATCH 12/25] Fixes in tests --- src/Algebra.jl | 6 ++++-- test/PLaplacianTests.jl | 7 ------- test/StokesOpenBoundaryTests.jl | 2 +- 3 files changed, 5 insertions(+), 10 deletions(-) diff --git a/src/Algebra.jl b/src/Algebra.jl index 47d6ca40..2f1bf1fd 100644 --- a/src/Algebra.jl +++ b/src/Algebra.jl @@ -169,9 +169,11 @@ function i_am_in(comm::MPIVoidVector) end function change_parts(x::Union{MPIArray,DebugArray,Nothing,MPIVoidVector}, new_parts; default=nothing) - x_new = map(new_parts) do _p - if isa(x,MPIArray) || isa(x,DebugArray) + x_new = map(new_parts) do p + if isa(x,MPIArray) PartitionedArrays.getany(x) + elseif isa(x,DebugArray) && (p <= length(x.items)) + x.items[p] else default end diff --git a/test/PLaplacianTests.jl b/test/PLaplacianTests.jl index 63a2d447..2347dfc4 100644 --- a/test/PLaplacianTests.jl +++ b/test/PLaplacianTests.jl @@ -7,13 +7,6 @@ using PartitionedArrays using Test using SparseArrays -# Overload required for the tests below -function Base.copy(a::PSparseMatrix) - a_matrix_partition = similar(a.matrix_partition) - copy!(a_matrix_partition, a.matrix_partition) - PSparseMatrix(a_matrix_partition,a.row_partition,a.col_partition) -end - function main(distribute,parts) main(distribute,parts,FullyAssembledRows(),SparseMatrixCSR{0,Float64,Int}) main(distribute,parts,SubAssembledRows(),SparseMatrixCSC{Float64,Int}) diff --git a/test/StokesOpenBoundaryTests.jl b/test/StokesOpenBoundaryTests.jl index d3a61b30..ca4213b5 100644 --- a/test/StokesOpenBoundaryTests.jl +++ b/test/StokesOpenBoundaryTests.jl @@ -3,7 +3,7 @@ module StokesOpenBoundaryTests using Gridap using LinearAlgebra using Test -using Gridap.ODEs +using Gridap.ODEs, Gridap.Algebra using GridapDistributed using PartitionedArrays From d01ed48280ccf5fcc41a5e17d6559a29f9be4fc7 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 22 Mar 2024 10:22:37 +1100 Subject: [PATCH 13/25] Updated compats to Gridap v0.18 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 511ee652..a883e0f0 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,7 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192" [compat] BlockArrays = "0.16.38" FillArrays = "0.8.4,1" -Gridap = "0.17.23" +Gridap = "0.18" LinearAlgebra = "1.3" MPI = "0.16, 0.17, 0.18, 0.19, 0.20" PartitionedArrays = "0.3.3" From b8f8457297d00e618f2ee70a194044fbc7698116 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 22 Mar 2024 10:31:12 +1100 Subject: [PATCH 14/25] Updated NEWS --- NEWS.md | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index c2a3145d..b67d6d61 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,13 +5,21 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.3.6] 2024-01-28 +## [0.4.0] 2024-03-22 + +### Changed + +- `DistributedCellField` now inherits from `CellField`. To accomodate the necessary API, we now save a pointer to the `DistributedTriangulation` where it is defined. This also requires `DistributedSingleFieldFESpace` to save the triangulation. Since PR[#141](https://github.com/gridap/GridapDistributed.jl/pull/141). +- All the distributed `Multifield` cellfield types are now represented by a `DistributedMultiFieldCellField`. Both `DistributedMultiFieldFEFunction` and `DistributedMultiFieldFEBasis` structs have been removed and replaced with constant aliases, which makes it more consistent with single-field types. Since PR[#141](https://github.com/gridap/GridapDistributed.jl/pull/141). +- Major refactor of ODE module. Implementation has been significantly simplified, while increasing the capability of the API. All `TransientDistributedObjects` structs have been removed, and replaced by `DistributedTransientObjects = DistributedObjects{TransientObject}`. Full support for EX/IM/IMEX methods. See Gridap's release for details. Since PR[#141](https://github.com/gridap/GridapDistributed.jl/pull/141). + +## [0.3.6] 2024-01-28 ### Added - Added redistribution for MultiFieldFESpaces. Since PR [#140](https://github.com/gridap/GridapDistributed.jl/pull/140). -### Fixed +### Fixed - Fixed issue [#142](https://github.com/gridap/GridapDistributed.jl/issues/142). Since PR [#142](https://github.com/gridap/GridapDistributed.jl/pull/142). From f1fee4169fdcb7490269cc50f2d3c7938586c175 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 22 Mar 2024 10:31:38 +1100 Subject: [PATCH 15/25] Bumped version to 0.4.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a883e0f0..9c0da0a4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridapDistributed" uuid = "f9701e48-63b3-45aa-9a63-9bc6c271f355" authors = ["S. Badia ", "A. F. Martin ", "F. Verdugo "] -version = "0.3.6" +version = "0.4.0" [deps] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" From 5f1e5c25000f6506fcc6a4c29a7a47bbeaf35bc0 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 27 Mar 2024 14:24:45 +1100 Subject: [PATCH 16/25] Added dot for MultiField --- src/Adaptivity.jl | 6 +++--- src/Geometry.jl | 4 ++++ src/MultiField.jl | 5 +++++ test/MultiFieldTests.jl | 8 ++++++++ 4 files changed, 20 insertions(+), 3 deletions(-) diff --git a/src/Adaptivity.jl b/src/Adaptivity.jl index 83313f86..d648260d 100644 --- a/src/Adaptivity.jl +++ b/src/Adaptivity.jl @@ -3,9 +3,9 @@ const DistributedAdaptedDiscreteModel{Dc,Dp} = GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:AdaptedDiscreteModel{Dc,Dp}}} -function DistributedAdaptedDiscreteModel(model ::DistributedDiscreteModel, - parent ::DistributedDiscreteModel, - glue ::AbstractArray{<:AdaptivityGlue}) +function DistributedAdaptedDiscreteModel(model :: DistributedDiscreteModel, + parent :: DistributedDiscreteModel, + glue :: AbstractArray{<:AdaptivityGlue}) models = map(local_views(model),local_views(parent),glue) do model, parent, glue AdaptedDiscreteModel(model,parent,glue) end diff --git a/src/Geometry.jl b/src/Geometry.jl index d80d9cb9..3a3b3568 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -398,6 +398,10 @@ function Geometry.get_background_model(a::DistributedTriangulation) a.model end +function Geometry.num_cells(a::DistributedTriangulation) + sum(map(trian->num_cells(trian),local_views(a))) +end + # Triangulation constructors function Geometry.Triangulation( diff --git a/src/MultiField.jl b/src/MultiField.jl index 56114b13..99bf8ad3 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -33,6 +33,11 @@ Base.iterate(m::DistributedMultiFieldCellField) = iterate(m.field_fe_fun) Base.iterate(m::DistributedMultiFieldCellField,state) = iterate(m.field_fe_fun,state) Base.getindex(m::DistributedMultiFieldCellField,field_id::Integer) = m.field_fe_fun[field_id] +function LinearAlgebra.dot(a::DistributedMultiFieldCellField,b::DistributedMultiFieldCellField) + @check num_fields(a) == num_fields(b) + return sum(map(dot,a.field_fe_fun,b.field_fe_fun)) +end + # DistributedMultiFieldFEFunction const DistributedMultiFieldFEFunction{A,B,T} = DistributedMultiFieldCellField{A,B,DistributedFEFunctionData{T}} diff --git a/test/MultiFieldTests.jl b/test/MultiFieldTests.jl index cb53e322..00f349ea 100644 --- a/test/MultiFieldTests.jl +++ b/test/MultiFieldTests.jl @@ -68,6 +68,14 @@ function main(distribute, parts, mfs) @test l2_error(p,ph1,dΩ) < 1.0e-9 @test l2_error(u,uh2,dΩ) < 1.0e-9 @test l2_error(p,ph2,dΩ) < 1.0e-9 + + a1(x,y) = ∫(x⋅y)dΩ + a2((u,p),(v,q)) = ∫(u⋅v + p⋅q)dΩ + A1 = assemble_matrix(a1,UxP,UxP) + A2 = assemble_matrix(a2,UxP,UxP) + + x = prandn(partition(axes(A1,2))) + @test norm(A1*x-A2*x) < 1.0e-9 end function main(distribute, parts) From e1e11f70d9405c6b9c6a524ec390bbc218b6fd6f Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 3 Apr 2024 14:51:35 +1100 Subject: [PATCH 17/25] Minor --- src/Algebra.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/Algebra.jl b/src/Algebra.jl index 2f1bf1fd..6a1a80f4 100644 --- a/src/Algebra.jl +++ b/src/Algebra.jl @@ -156,17 +156,18 @@ function get_part_id(comm::MPI.Comm) id end +""" + i_am_in(comm::MPIArray) + i_am_in(comm::DebugArray) + + Returns `true` if the processor is part of the subcommunicator `comm`. +""" function i_am_in(comm::MPI.Comm) get_part_id(comm) >=0 end - -function i_am_in(comm::MPIArray) - i_am_in(comm.comm) -end - -function i_am_in(comm::MPIVoidVector) - i_am_in(comm.comm) -end +@inline i_am_in(comm::MPIArray) = i_am_in(comm.comm) +@inline i_am_in(comm::MPIVoidVector) = i_am_in(comm.comm) +@inline i_am_in(comm::DebugArray) = true function change_parts(x::Union{MPIArray,DebugArray,Nothing,MPIVoidVector}, new_parts; default=nothing) x_new = map(new_parts) do p From 4fb1b1fbf2b7b986bca75ee6712b6cb4e4440529 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 3 Apr 2024 15:39:48 +1100 Subject: [PATCH 18/25] Minor --- src/Algebra.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/Algebra.jl b/src/Algebra.jl index 6a1a80f4..4af8c5c0 100644 --- a/src/Algebra.jl +++ b/src/Algebra.jl @@ -146,6 +146,9 @@ function num_parts(comm::MPI.Comm) end nparts end +@inline num_parts(comm::MPIArray) = num_parts(comm.comm) +@inline num_parts(comm::DebugArray) = length(comm.items) +@inline num_parts(comm::MPIVoidVector) = num_parts(comm.comm) function get_part_id(comm::MPI.Comm) if comm != MPI.COMM_NULL @@ -155,6 +158,8 @@ function get_part_id(comm::MPI.Comm) end id end +@inline get_part_id(comm::MPIArray) = get_part_id(comm.comm) +@inline get_part_id(comm::MPIVoidVector) = get_part_id(comm.comm) """ i_am_in(comm::MPIArray) From 73b87b2b987bb74c8a5b9c76b0ebb38b1bdbab15 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 9 Apr 2024 12:39:25 +1000 Subject: [PATCH 19/25] Minor --- src/MultiField.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/MultiField.jl b/src/MultiField.jl index 99bf8ad3..1a945bcc 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -113,6 +113,10 @@ function FESpaces.zero_dirichlet_values(f::DistributedMultiFieldFESpace) map(zero_dirichlet_values,f.field_fe_space) end +function FESpaces.get_dirichlet_dof_values(f::DistributedMultiFieldFESpace) + return map(get_dirichlet_dof_values,f.field_fe_space) +end + function FESpaces.FEFunction(f::DistributedMultiFieldFESpace,x::AbstractVector,isconsistent=false) free_values = change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true) part_fe_fun = map(FEFunction,f.part_fe_space,partition(free_values)) @@ -129,9 +133,10 @@ end function FESpaces.FEFunction( f::DistributedMultiFieldFESpace,x::AbstractVector, dirichlet_values::AbstractArray{<:AbstractVector},isconsistent=false - ) +) free_values = GridapDistributed.change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true) - part_fe_fun = map(FEFunction,f.part_fe_space,partition(free_values)) + part_dirvals = to_parray_of_arrays(dirichlet_values) + part_fe_fun = map(FEFunction,f.part_fe_space,partition(free_values),part_dirvals) field_fe_fun = DistributedSingleFieldFEFunction[] for i in 1:num_fields(f) free_values_i = restrict_to_field(f,free_values,i) From e2ecabfb585c319b5743317c2b789bed35f3c1df Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 12 Apr 2024 12:26:13 +1000 Subject: [PATCH 20/25] Bugfix: Block indexing allocates copies --- NEWS.md | 2 +- src/Algebra.jl | 2 +- src/BlockPartitionedArrays.jl | 11 +++++++---- src/MultiField.jl | 5 +++-- 4 files changed, 12 insertions(+), 8 deletions(-) diff --git a/NEWS.md b/NEWS.md index b67d6d61..b109d014 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,7 +5,7 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). -## [0.4.0] 2024-03-22 +## [0.4.0] 2024-04-12 ### Changed diff --git a/src/Algebra.jl b/src/Algebra.jl index 4af8c5c0..cd1a1634 100644 --- a/src/Algebra.jl +++ b/src/Algebra.jl @@ -373,7 +373,7 @@ end function local_views(a::BlockPMatrix,new_rows::BlockPRange,new_cols::BlockPRange) vals = map(CartesianIndices(blocksize(a))) do I - local_views(a[Block(I)],new_rows[Block(I[1])],new_cols[Block(I[2])]) + local_views(blocks(a)[I],blocks(new_rows)[I],blocks(new_cols)[I]) end |> to_parray_of_arrays return map(mortar,vals) end diff --git a/src/BlockPartitionedArrays.jl b/src/BlockPartitionedArrays.jl index 5661432f..cd9e1c6e 100644 --- a/src/BlockPartitionedArrays.jl +++ b/src/BlockPartitionedArrays.jl @@ -341,16 +341,19 @@ function Base.broadcasted(f, a::Union{BlockPArray,BlockPBroadcasted}, b::Number) return BlockPBroadcasted(blocks_out,blockaxes(a)) end -function Base.broadcasted(f, - a::Union{BlockPArray,BlockPBroadcasted}, - b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}) +function Base.broadcasted( + f, + a::Union{BlockPArray,BlockPBroadcasted}, + b::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}} +) Base.broadcasted(f,a,Base.materialize(b)) end function Base.broadcasted( f, a::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}}, - b::Union{BlockPArray,BlockPBroadcasted}) + b::Union{BlockPArray,BlockPBroadcasted} +) Base.broadcasted(f,Base.materialize(a),b) end diff --git a/src/MultiField.jl b/src/MultiField.jl index 1a945bcc..e8927d92 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -101,8 +101,9 @@ function FESpaces.get_free_dof_ids(fs::DistributedMultiFieldFESpace) end function MultiField.restrict_to_field( - f::DistributedMultiFieldFESpace,free_values::AbstractVector,field::Integer) - values = map(f.part_fe_space,partition(free_values)) do u,x + f::DistributedMultiFieldFESpace,free_values::AbstractVector,field::Integer +) + values = map(local_views(f),partition(free_values)) do u,x restrict_to_field(u,x,field) end gids = f.field_fe_space[field].gids From a3430324cda3981dd6d14253ba9e6979677be9a6 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 12 Apr 2024 17:45:10 +1000 Subject: [PATCH 21/25] Fixed BlockPArray bugs --- src/BlockPartitionedArrays.jl | 54 +++++++++++++++++++++++++++++------ test/MultiFieldTests.jl | 8 +++--- 2 files changed, 49 insertions(+), 13 deletions(-) diff --git a/src/BlockPartitionedArrays.jl b/src/BlockPartitionedArrays.jl index cd9e1c6e..f173e5ce 100644 --- a/src/BlockPartitionedArrays.jl +++ b/src/BlockPartitionedArrays.jl @@ -145,17 +145,19 @@ end function Base.copyto!(y::BlockPVector,x::BlockPVector) @check blocklength(x) == blocklength(y) - for i in blockaxes(x,1) - copyto!(y[i],x[i]) + yb, xb = blocks(y), blocks(x) + for i in 1:blocksize(x,1) + copyto!(yb[i],xb[i]) end return y end function Base.copyto!(y::BlockPMatrix,x::BlockPMatrix) @check blocksize(x) == blocksize(y) - for i in blockaxes(x,1) - for j in blockaxes(x,2) - copyto!(y[i,j],x[i,j]) + yb, xb = blocks(y), blocks(x) + for i in 1:blocksize(x,1) + for j in 1:blocksize(x,2) + copyto!(yb[i,j],xb[i,j]) end end return y @@ -169,6 +171,8 @@ function Base.fill!(a::BlockPVector,v) end function Base.sum(a::BlockPArray) + # TODO: This could use a single communication, instead of one for each block + # TODO: We could implement a generic reduce, that we apply to sum, all, any, etc.. return sum(map(sum,blocks(a))) end @@ -284,15 +288,47 @@ end # LinearAlgebra API +function Base.:*(a::Number,b::BlockArray) + mortar(map(bi -> a*bi,blocks(b))) +end +Base.:*(b::BlockPMatrix,a::Number) = a*b +Base.:/(b::BlockPVector,a::Number) = (1/a)*b + +function Base.:*(a::BlockPMatrix,b::BlockPVector) + c = similar(b) + mul!(c,a,b) + return c +end + +for op in (:+,:-) + @eval begin + function Base.$op(a::BlockPArray) + mortar(map($op,blocks(a))) + end + function Base.$op(a::BlockPArray,b::BlockPArray) + @assert blocksize(a) == blocksize(b) + mortar(map($op,blocks(a),blocks(b))) + end + end +end + function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector) + o = one(eltype(A)) + mul!(y,A,x,o,o) +end + +function LinearAlgebra.mul!(y::BlockPVector,A::BlockPMatrix,x::BlockPVector,α::Number,β::Number) + yb, Ab, xb = blocks(y), blocks(A), blocks(x) z = zero(eltype(y)) o = one(eltype(A)) - for i in blockaxes(A,2) - fill!(y[i],z) - for j in blockaxes(A,2) - mul!(y[i],A[i,j],x[j],o,o) + for i in 1:blocksize(A,1) + fill!(yb[i],z) + for j in 1:blocksize(A,2) + mul!(yb[i],Ab[i,j],xb[j],α,o) end + rmul!(yb[i],β) end + return y end function LinearAlgebra.dot(x::BlockPVector,y::BlockPVector) diff --git a/test/MultiFieldTests.jl b/test/MultiFieldTests.jl index 00f349ea..9a0475ab 100644 --- a/test/MultiFieldTests.jl +++ b/test/MultiFieldTests.jl @@ -1,8 +1,7 @@ module MultiFieldTests using Gridap -using Gridap.FESpaces -using Gridap.MultiField +using Gridap.FESpaces, Gridap.MultiField, Gridap.Algebra using GridapDistributed using PartitionedArrays using Test @@ -74,8 +73,9 @@ function main(distribute, parts, mfs) A1 = assemble_matrix(a1,UxP,UxP) A2 = assemble_matrix(a2,UxP,UxP) - x = prandn(partition(axes(A1,2))) - @test norm(A1*x-A2*x) < 1.0e-9 + x1 = allocate_in_domain(A1); fill!(x1,1.0) + x2 = allocate_in_domain(A2); fill!(x2,1.0) + @test norm(A1*x1-A2*x2) < 1.0e-9 end function main(distribute, parts) From 3d56e8ab8bd177695c6d6e47978e60a41af41e55 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 12 Apr 2024 17:56:17 +1000 Subject: [PATCH 22/25] Changed sign of residual in tests --- test/HeatEquationTests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/HeatEquationTests.jl b/test/HeatEquationTests.jl index 6bda0178..063cb24f 100644 --- a/test/HeatEquationTests.jl +++ b/test/HeatEquationTests.jl @@ -38,7 +38,7 @@ function main(distribute,parts) assembler = SparseMatrixAssembler(U,V0,SubAssembledRows()) op_constant = TransientLinearFEOperator( - (a,m),(t,v) -> (-1)*b(t,v),U,V0,constant_forms=(true,true),assembler=assembler + (a,m),b,U,V0,constant_forms=(true,true),assembler=assembler ) t0 = 0.0 From d6d52001b6f46941e5c5ee5d325d6f50273fd315 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sun, 14 Apr 2024 17:30:44 +1000 Subject: [PATCH 23/25] Minor --- test/HeatEquationTests.jl | 7 ++++--- test/StokesOpenBoundaryTests.jl | 23 +++++++++++++---------- 2 files changed, 17 insertions(+), 13 deletions(-) diff --git a/test/HeatEquationTests.jl b/test/HeatEquationTests.jl index 063cb24f..e86f0d9d 100644 --- a/test/HeatEquationTests.jl +++ b/test/HeatEquationTests.jl @@ -10,9 +10,10 @@ function main(distribute,parts) ranks = distribute(LinearIndices((prod(parts),))) θ = 0.2 - u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t - u(t::Real) = x -> u(x,t) - f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x) + ut(t) = x -> (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t + u = TimeSpaceFunction(ut) + ft(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x) + f = TimeSpaceFunction(ft) domain = (0,1,0,1) model = CartesianDiscreteModel(ranks,parts,domain,(4,4)) diff --git a/test/StokesOpenBoundaryTests.jl b/test/StokesOpenBoundaryTests.jl index ca4213b5..a2647645 100644 --- a/test/StokesOpenBoundaryTests.jl +++ b/test/StokesOpenBoundaryTests.jl @@ -12,16 +12,19 @@ function main(distribute,parts) θ = 0.5 - u(x,t) = VectorValue(x[1],x[2])*t - u(t::Real) = x -> u(x,t) - - p(x,t) = (x[1]-x[2])*t - p(t::Real) = x -> p(x,t) - q(x) = t -> p(x,t) - - f(t) = x -> ∂t(u)(t)(x)-Δ(u(t))(x)+ ∇(p(t))(x) - g(t) = x -> (∇⋅u(t))(x) - h(t) = x -> ∇(u(t))(x)⋅VectorValue(0.0,1.0) - p(t)(x)*VectorValue(0.0,1.0) + ut(t) = x -> VectorValue(x[1],x[2])*t + u = TimeSpaceFunction(ut) + + pt(t) = x -> (x[1]-x[2])*t + p = TimeSpaceFunction(pt) + q(x) = t -> p(t,x) + + ft(t) = x -> ∂t(u)(t)(x)-Δ(u(t))(x)+ ∇(p(t))(x) + gt(t) = x -> (∇⋅u(t))(x) + ht(t) = x -> ∇(u(t))(x)⋅VectorValue(0.0,1.0) - p(t)(x)*VectorValue(0.0,1.0) + f = TimeSpaceFunction(ft) + g = TimeSpaceFunction(gt) + h = TimeSpaceFunction(ht) domain = (0,1,0,1) partition = (4,4) From 1082ae6e945f8280c47048d49802839333da1fb2 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 15 Apr 2024 08:53:18 +1000 Subject: [PATCH 24/25] Fixed tests after 0.18.1 --- test/HeatEquationTests.jl | 2 +- test/StokesOpenBoundaryTests.jl | 14 +++++--------- 2 files changed, 6 insertions(+), 10 deletions(-) diff --git a/test/HeatEquationTests.jl b/test/HeatEquationTests.jl index e86f0d9d..fa96dc11 100644 --- a/test/HeatEquationTests.jl +++ b/test/HeatEquationTests.jl @@ -12,7 +12,7 @@ function main(distribute,parts) ut(t) = x -> (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*t u = TimeSpaceFunction(ut) - ft(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x) + ft(t) = x -> ∂t(u)(t,x) - Δ(u)(t,x) f = TimeSpaceFunction(ft) domain = (0,1,0,1) diff --git a/test/StokesOpenBoundaryTests.jl b/test/StokesOpenBoundaryTests.jl index a2647645..3af86558 100644 --- a/test/StokesOpenBoundaryTests.jl +++ b/test/StokesOpenBoundaryTests.jl @@ -19,9 +19,9 @@ function main(distribute,parts) p = TimeSpaceFunction(pt) q(x) = t -> p(t,x) - ft(t) = x -> ∂t(u)(t)(x)-Δ(u(t))(x)+ ∇(p(t))(x) - gt(t) = x -> (∇⋅u(t))(x) - ht(t) = x -> ∇(u(t))(x)⋅VectorValue(0.0,1.0) - p(t)(x)*VectorValue(0.0,1.0) + ft(t) = x -> ∂t(u)(t,x) - Δ(u)(t,x) + ∇(p)(t,x) + gt(t) = x -> (∇⋅u)(t,x) + ht(t) = x -> ∇(u)(t,x)⋅VectorValue(0.0,1.0) - p(t,x)*VectorValue(0.0,1.0) f = TimeSpaceFunction(ft) g = TimeSpaceFunction(gt) h = TimeSpaceFunction(ht) @@ -99,15 +99,11 @@ function main(distribute,parts) ph_tn = xh_tn[2] #writevtk(Ω,"output/tmp_stokes_OB_sol_$tn.vtu",cellfields=["u"=>uh_tn,"p"=>ph_tn]) e = u(tn) - uh_tn - el2 = sqrt(sum( ∫(l2(e))dΩ )) + el2 = sqrt(sum(∫(l2(e))dΩ)) e = p(tn) - ph_tn - el2 = sqrt(sum( ∫(l2(e))dΩ )) + el2 = sqrt(sum(∫(l2(e))dΩ)) @test el2 < tol end end -with_debug() do distribute - main(distribute,(2,2)) -end - end #module From 6a5a51f49151ef773fe6ec6bdf76c38f2355f039 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 15 Apr 2024 09:09:17 +1000 Subject: [PATCH 25/25] Small typo --- src/BlockPartitionedArrays.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/BlockPartitionedArrays.jl b/src/BlockPartitionedArrays.jl index f173e5ce..25e43466 100644 --- a/src/BlockPartitionedArrays.jl +++ b/src/BlockPartitionedArrays.jl @@ -288,7 +288,7 @@ end # LinearAlgebra API -function Base.:*(a::Number,b::BlockArray) +function Base.:*(a::Number,b::BlockPArray) mortar(map(bi -> a*bi,blocks(b))) end Base.:*(b::BlockPMatrix,a::Number) = a*b