From 139f6fcbc333910dde3c3e2f6a3daea69fdfe0ef Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 29 Jul 2024 16:26:19 +1000 Subject: [PATCH 1/6] Fixed interpolators for VectorValues --- src/CellData.jl | 111 +++++++++++++++++++++++++----------------- test/CellDataTests.jl | 26 ++++++---- 2 files changed, 82 insertions(+), 55 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index 95518ce..6aded00 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -358,70 +358,91 @@ end Arrays.evaluate!(cache,f::DistributedCellField,x::Point) = evaluate!(cache,Interpolable(f),x) Arrays.evaluate!(cache,f::DistributedCellField,x::AbstractVector{<:Point}) = evaluate!(cache,Interpolable(f),x) -struct DistributedInterpolable{A} <: Function +struct DistributedInterpolable{Tx,Ty,A} <: Function interps::A + function DistributedInterpolable(interps::AbstractArray{<:Interpolable}) + Tx,Ty = map(interps) do I + fi = I.uh + trian = get_triangulation(fi) + x = mean(testitem(get_cell_coordinates(trian))) + return typeof(x), return_type(fi,x) + end |> tuple_of_arrays + Tx = getany(Tx) + Ty = getany(Ty) + A = typeof(interps) + new{Tx,Ty,A}(interps) + end end + local_views(a::DistributedInterpolable) = a.interps function Interpolable(f::DistributedCellField;kwargs...) - interps = map(local_views(f)) do fun - Interpolable(fun,kwargs...) + interps = map(local_views(f)) do f + Interpolable(f,kwargs...) end DistributedInterpolable(interps) end (a::DistributedInterpolable)(x) = evaluate(a,x) -function Gridap.Arrays.return_cache(I::DistributedInterpolable,x::Point) - caches = map(local_views(I)) do fi - trian = get_triangulation(fi.uh) - y=mean(testitem(get_cell_coordinates(trian))) - @check typeof(testitem(x)) == typeof(y) "Can only evaluate DistributedInterpolable at physical points of the same dimension of the underlying triangulation" - return_cache(fi,y) +Arrays.return_cache(f::DistributedInterpolable,x::Point) = Arrays.return_cache(f,[x]) + +Arrays.evaluate!(caches,I::DistributedInterpolable,x::Point) = evaluate!(caches,I,[x])[1] + +function Arrays.return_cache(I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty} + msg = "Can only evaluate DistributedInterpolable at physical points of the same dimension of the underlying triangulation" + @check Tx == eltype(x) msg + caches = map(local_views(I)) do I + trian = get_triangulation(I.uh) + y = mean(testitem(get_cell_coordinates(trian))) + return_cache(I,y) + end + caches +end + +function Gridap.Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty} + infty(::Type{T}) where T = -T(Inf) + infty(::Type{VectorValue{D,T}}) where {D,T} = VectorValue(fill(infty(T),D)...) + combine(a,b) = map(max,a,b) + function array_reduce(f,a) + b = gather(a) + c = map_main(b) do b + c = fill(infty(eltype((first(b)))),length(first(b))) + for bi in b + c = f(c,bi) + end + return c + end + return getany(emit(c)) end - caches -end -Gridap.Arrays.return_cache(f::DistributedInterpolable,x::AbstractVector{<:Point}) = Gridap.Arrays.return_cache(f,testitem(x)) -function Gridap.Arrays.evaluate!(caches,I::DistributedInterpolable,x::Point) - y=map(local_views(I),local_views(caches)) do fi,cache + # Evaluate in local portions of the domain. Set to -Inf if the point is not found. + nx = length(x) + y = map(local_views(I),local_views(cache)) do I, cache + y = Vector{Ty}(undef,nx) + for (i,xi) in enumerate(x) try - evaluate!(cache,fi,x) + y[i] = evaluate!(cache,I,xi) catch - -Inf + y[i] = infty(Ty) end - end - # reduce(max,y) - z=gather(y) - map_main(local_views(z)) do zi - reduce(max,zi) end -end + return y + end -function Gridap.Arrays.evaluate!(caches,I::DistributedInterpolable,v::AbstractVector{<:Point}) - n=length(local_views(I)) - m=length(v) - y=map(local_views(I),local_views(caches)) do fi,cache - w=Vector{Float64}(undef,m) - for (i,x) in enumerate(v) - try - w[i]=evaluate!(cache,fi,x) - catch - w[i]=-Inf - end - end - return w - end - # z=gather(y,destination=:all) - z=gather(y) - map_main(local_views(z)) do zi - w=Vector{Float64}(undef,m) - for i=0:m-1 - w[i+1]=reduce(max,zi.data[zi.ptrs[1:n].+i]) - end - return w + # Combine the results + if Ty <: VectorValue + w_d = Vector{Vector{eltype(Ty)}}(undef,num_components(Ty)) + for d in 1:num_components(Ty) + y_d = map(y_p -> map(y_p_i -> y_p_i[d],y_p),y) + w_d[d] = array_reduce(combine,y_d) end - # reduce((v,w)->broadcast(max,v,w),y) + w = map(VectorValue,w_d...) + else + w = array_reduce(combine,y) + end + + return w end # Support for distributed Dirac deltas diff --git a/test/CellDataTests.jl b/test/CellDataTests.jl index 2485374..3baa2b8 100644 --- a/test/CellDataTests.jl +++ b/test/CellDataTests.jl @@ -58,7 +58,7 @@ function main(distribute,parts) u3 = CellField(2.0,Ω) u = _my_op∘(u1,u2,u3) - order = 1 + order = 1 reffe = ReferenceFE(lagrangian,Float64,order) V = TestFESpace(model,reffe) uh = interpolate_everywhere(x->x[1]+x[2],V) @@ -67,15 +67,21 @@ function main(distribute,parts) x3 = Point(0.9,0.9) v = [x1,x2,x3] - u1 = uh(x1) - u2 = uh(x2) - uv = uh(v) + @test uh(x1) == 0.2 + @test uh(x2) == 1.0 + @test uh(v) == [0.2,1.0,1.8] - map_main(u1,u2,uv) do u1,u2,v - @test u1 == 0.2 - @test u2 == 1.0 - @test v ==[0.2,1.0,1.8] - end + reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) + V = TestFESpace(model,reffe) + uh = interpolate_everywhere(x->x,V) + x1 = Point(0.1,0.1) + x2 = Point(0.1,0.9) + x3 = Point(0.9,0.9) + v = [x1,x2,x3] + + @test uh(x1) == x1 + @test uh(x2) == x2 + @test uh(v) == v # Point δ δ=DiracDelta{0}(model;tags=2) @@ -89,7 +95,7 @@ function main(distribute,parts) @test sum(δ(f)) ≈ 8.0 @test sum(δ(3.0)) ≈ 12.0 @test sum(δ(x->2*x)) ≈ VectorValue(16.0,0.0) - + end end # module From ac12ab3dec2a20466e252f5618091fa8b27a29a7 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 29 Jul 2024 16:29:34 +1000 Subject: [PATCH 2/6] Updated NEWS.md --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index 24ebcf1..65ba24d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,12 @@ 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). +## Unreleased + +### Fixed + +- Fixed distributed interpolators for Vector-Valued FESpaces. Since PR[#152](https://github.com/gridap/GridapDistributed.jl/pull/152). + ## [0.4.3] 2024-07-18 ### Added From d8ab3b69c6aee702dda5a43df65db558d3970ccc Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Mon, 29 Jul 2024 23:08:32 +1000 Subject: [PATCH 3/6] Minor --- src/CellData.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index 6aded00..6994ea5 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -355,8 +355,8 @@ f(s) instead of s(f)? end # Interpolation at arbitrary points (returns -Inf if the point is not found) -Arrays.evaluate!(cache,f::DistributedCellField,x::Point) = evaluate!(cache,Interpolable(f),x) -Arrays.evaluate!(cache,f::DistributedCellField,x::AbstractVector{<:Point}) = evaluate!(cache,Interpolable(f),x) +Arrays.evaluate!(cache,f::DistributedCellField,x::Point) = evaluate(Interpolable(f),x) +Arrays.evaluate!(cache,f::DistributedCellField,x::AbstractVector{<:Point}) = evaluate(Interpolable(f),x) struct DistributedInterpolable{Tx,Ty,A} <: Function interps::A @@ -385,9 +385,8 @@ end (a::DistributedInterpolable)(x) = evaluate(a,x) -Arrays.return_cache(f::DistributedInterpolable,x::Point) = Arrays.return_cache(f,[x]) - -Arrays.evaluate!(caches,I::DistributedInterpolable,x::Point) = evaluate!(caches,I,[x])[1] +Arrays.return_cache(f::DistributedInterpolable,x::Point) = return_cache(f,[x]) +Arrays.evaluate!(caches,I::DistributedInterpolable,x::Point) = first(evaluate!(caches,I,[x])) function Arrays.return_cache(I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty} msg = "Can only evaluate DistributedInterpolable at physical points of the same dimension of the underlying triangulation" @@ -400,7 +399,7 @@ function Arrays.return_cache(I::DistributedInterpolable{Tx,Ty},x::AbstractVector caches end -function Gridap.Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty} +function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty} infty(::Type{T}) where T = -T(Inf) infty(::Type{VectorValue{D,T}}) where {D,T} = VectorValue(fill(infty(T),D)...) combine(a,b) = map(max,a,b) From a00eaa49c8f750e3e93ee227e561083c62bfefdb Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 30 Jul 2024 13:56:20 +1000 Subject: [PATCH 4/6] Fix --- src/CellData.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index 6994ea5..9ceb015 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -404,9 +404,10 @@ function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVec infty(::Type{VectorValue{D,T}}) where {D,T} = VectorValue(fill(infty(T),D)...) combine(a,b) = map(max,a,b) function array_reduce(f,a) + T = getany(map(eltype,a)) b = gather(a) - c = map_main(b) do b - c = fill(infty(eltype((first(b)))),length(first(b))) + c = map_main(b;otherwise=b->T[]) do b + c = fill(infty(T),length(first(b))) for bi in b c = f(c,bi) end From e2ec4de6ef81e4ac5498859bc50781618fc7871c Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 30 Jul 2024 14:36:20 +1000 Subject: [PATCH 5/6] Reduce memory and communication footprint --- src/CellData.jl | 61 ++++++++++++++++++++++++++----------------------- 1 file changed, 33 insertions(+), 28 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index 9ceb015..36f0787 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -400,46 +400,51 @@ function Arrays.return_cache(I::DistributedInterpolable{Tx,Ty},x::AbstractVector end function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty} - infty(::Type{T}) where T = -T(Inf) - infty(::Type{VectorValue{D,T}}) where {D,T} = VectorValue(fill(infty(T),D)...) - combine(a,b) = map(max,a,b) - function array_reduce(f,a) - T = getany(map(eltype,a)) - b = gather(a) - c = map_main(b;otherwise=b->T[]) do b - c = fill(infty(T),length(first(b))) - for bi in b - c = f(c,bi) - end - return c - end - return getany(emit(c)) - end + _allgather(x) = getdata(getany(gather(x;destination=:all))) - # Evaluate in local portions of the domain. Set to -Inf if the point is not found. + # Evaluate in local portions of the domain. Only keep points inside the domain. nx = length(x) - y = map(local_views(I),local_views(cache)) do I, cache - y = Vector{Ty}(undef,nx) + my_ids, my_vals = map(local_views(I),local_views(cache)) do I, cache + ids = Vector{Int}(undef,nx) + vals = Vector{Ty}(undef,nx) + k = 1 + yi = zero(Ty) for (i,xi) in enumerate(x) + inside = true try - y[i] = evaluate!(cache,I,xi) + yi = evaluate!(cache,I,xi) catch - y[i] = infty(Ty) + inside = false + end + if inside + ids[k] = i + vals[k] = yi + k += 1 end end - return y + resize!(ids,k-1) + resize!(vals,k-1) + return ids, vals end - # Combine the results + # Communicate results, so that every (id,value) pair is known by every process if Ty <: VectorValue - w_d = Vector{Vector{eltype(Ty)}}(undef,num_components(Ty)) - for d in 1:num_components(Ty) - y_d = map(y_p -> map(y_p_i -> y_p_i[d],y_p),y) - w_d[d] = array_reduce(combine,y_d) + D = num_components(Ty) + vals_d = Vector{Vector{eltype(Ty)}}(undef,D) + for d in 1:D + my_vals_d = map(y_p -> map(y_p_i -> y_p_i[d],y_p),my_vals) + vals_d[d] = _allgather(my_vals_d) end - w = map(VectorValue,w_d...) + vals = map(VectorValue,w_d...) else - w = array_reduce(combine,y) + vals = _allgather(my_vals) + end + ids = _allgather(my_ids) + + # Combine results + w = Vector{Ty}(undef,nx) + for (i,v) in zip(ids,vals) + w[i] = v end return w From d7ea641df39e193dec535c3b4ac443be56bbb4be Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 30 Jul 2024 15:40:56 +1000 Subject: [PATCH 6/6] Fixes --- src/CellData.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/CellData.jl b/src/CellData.jl index 36f0787..6befc88 100644 --- a/src/CellData.jl +++ b/src/CellData.jl @@ -400,7 +400,7 @@ function Arrays.return_cache(I::DistributedInterpolable{Tx,Ty},x::AbstractVector end function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVector{<:Point}) where {Tx,Ty} - _allgather(x) = getdata(getany(gather(x;destination=:all))) + _allgather(x) = PartitionedArrays.getdata(getany(gather(x;destination=:all))) # Evaluate in local portions of the domain. Only keep points inside the domain. nx = length(x) @@ -425,7 +425,7 @@ function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVec resize!(ids,k-1) resize!(vals,k-1) return ids, vals - end + end |> tuple_of_arrays # Communicate results, so that every (id,value) pair is known by every process if Ty <: VectorValue @@ -435,7 +435,7 @@ function Arrays.evaluate!(cache,I::DistributedInterpolable{Tx,Ty},x::AbstractVec my_vals_d = map(y_p -> map(y_p_i -> y_p_i[d],y_p),my_vals) vals_d[d] = _allgather(my_vals_d) end - vals = map(VectorValue,w_d...) + vals = map(VectorValue,vals_d...) else vals = _allgather(my_vals) end