Interpolate everywhere arbitrary points #632
Conversation
…erpolate_everywhere_arbitrary_points
…erpolate_everywhere_arbitrary_points
Codecov Report
@@ Coverage Diff @@
## master #632 +/- ##
==========================================
+ Coverage 87.66% 87.70% +0.04%
==========================================
Files 134 134
Lines 14267 14304 +37
==========================================
+ Hits 12507 12546 +39
+ Misses 1760 1758 -2
Continue to review full report at Codecov.
|
|
Hi @Balaje ! I'll review the PR later in more depth. At first sight, I have a couple of questions:
|
|
Hi, @fverdugo,
We want to interpolate a FE function from a FE space OLD to a FE space NEW, OLD and NEW defined in different triangulations. At first, we thought it could be a good idea to have fₕ = interpolate_everywhere(fe_function_in_OLD_fe_space,NEW_fe_space)But fe_function_in_NEW_fe_space = CellField(NEW_fe_space,fe_function_in_OLD_fe_space) # New method!
fₕ = interpolate_everywhere(fe_function_in_NEW_fe_space,NEW_fe_space)This solution does not touch the
To create the FE function in the We just need to figure out the best way to implement the following functionality: Given Maybe we could start from here. |
|
The key part in the interpolation routines is this helper function Gridap.jl/src/FESpaces/SingleFieldFESpaces.jl Line 303 in cdb45ce function _cell_vals(fs::SingleFieldFESpace,object)
s = get_fe_dof_basis(fs)
trian = get_triangulation(s)
f = CellField(object,trian,DomainStyle(s))
cell_vals = s(f)
endIn your case, Gridap.jl/src/CellData/CellFields.jl Line 95 in cdb45ce function CellField(f::CellField,trian::Triangulation,domain_style::DomainStyle)
change_domain(f,trian,domain_style)
endThus, the new case to be implemented is the call to Fixing I would start by refactoring you code so that we rely on |
|
In a second stage, we can discuss about the signature |
|
Hi @fverdugo,
This was precisely @Balaje's first approach. Which, indeed, simplifies quite a bit the implementation. However, this is dangerous since we would effectively loose the following error when accidentally using incompatible triangulations: If I'm not mistaken, then, if I use a triangulation that is not compatible, it would always interpolate instead of throwing the error. Is that something that we want? |
|
@oriolcg good question! I see pros and cos. It would be more "flexible" if the triangulations are not compatible try to interpolate as a fall-back, but it is also more dangerous since it is more difficult to guarantee that you are not interpolating.... Let me think about it. |
|
What about this? # User API
#
# by default, do not allow interpolation (just like now)
uh_new = interpolate_everywhere(uh_old,Vh_new)
# If you want to allow interpolation, pass this as an option
options = InterpolationOptions(interpolate=true)
uh_new = interpolate_everywhere(uh_old,Vh_new,options)
# we could even add more params into the options struct.
# It is likely that you need them at some point
# e.g. for passing tolerances
options = InterpolationOptions(
interpolate=true,
search_method=:kdtree,
tol=1e-6)
uh_new = interpolate_everywhere(uh_old,Vh_new,options)
# we could follow the same pattern for `interpolate` and `interpolate_dirichlet`.
# it also make sense for them.
# Possible implementation (incomplete)
struct InterpolationOptions
interpolation::Bool
search_method::Symbol
tol::Float64
function InterpolationOptions(;
interpolation=true,search_method=:kdtee,tol=1e-6)
new(interpolation,search_method,tol)
end
end
# We can assume that f::SingleFieldFEFunction
function interpolate_everywhere(
f::SingleFieldFEFunction,V::SingleFieldFESpace,options::InterpolationOptions)
if !options.interpolate
interpolate_everywhere(f,V)
else
s = get_fe_dof_basis(V)
trian = get_triangulation(s)
f_new = change_domain(f,trian,DomainStyle(s),options) # this is the new method to implement
cell_vals = s(f_new)
end
end
function change_domain(
f::SingleFieldFEFunction,trian::Triangulation,ds::DomainStyle,options::InterpolationOptions)
if !options.interpolate
change_domain(f,trian,ds)
else
# your interpolation stuff here
end
end
Perhaps we can find a shorter name for |
|
Hi, @fverdugo, upon reviewing the refactoring you propose, I realised that we can do: interpolate_everywhere(x -> uh_old(x),Vh_new)without touching Now, the question is: uh_new ?= interpolate_everywhere(x -> uh_old(x),Vh_new)I think it is, at least, @Balaje tests pass if I use Given this, I wonder if we can further simplify the refactor above or I am overlooking something... |
@ericneiva this should work, but I am afraid that it would be very slow. If I am not wrong, this is equivalent to generate the kde tree at each gauss point. |
In any case, I find this very elegant. Perhaps there is a way to speed up it and also to pass some options for the searching routines. Something like interpolate_everywhere(Interpolable(uh_old),Vh_new)
# or
interpolate_everywhere(Interpolable(uh_old,tol=1e-6),Vh_new)I like more this last option than adding an extra argument to |
|
We are closer. This looks quite good IMO: # by defining this
struct Interpolable <: Function
uh::A
tol::Float64
function Interpolable(uh;tol=1e-6)
new{typeof(uh)}(uh,tol)
end
end
(f::Interpolable)(x) = f.uh(x)
# this should work out of the box (but very inefficienlty)
uh_new = interpolate(Interpolable(uh_old),Vh_new)
uh_new = interpolate_everywhere(Interpolable(uh_old),Vh_new)
uh_new = interpolate_dirichlet(Interpolable(uh_old),Vh_new)
# To speed up things we can improve _cell_vals
function _cell_vals(V::SingleFieldFESpace,f::Interpolable)
# Interpolation stuff here.
# It can use parameters inside stored inside f.
# The result should be the cell-wise dof values
# in the new space.
# Gridap will handle the rest.
end
|
|
Hi @fverdugo I had a chat with @ericneiva and @santiagobadia today, and we agreed to stick with the implementation that involves transforming the I have the implementation separately here - /~https://github.com/Balaje/GSoC-2021/blob/main/Interpolation/interpolable.jl, and I'll push the changes to my local repository shortly. Thanks. |
|
Yes, moving the dofs to the physical space is another option, which for Lagrangian and moment-based dof bases can be done without using the inverse map. A couple of comments regarding /~https://github.com/Balaje/GSoC-2021/blob/main/Interpolation/interpolable.jl This implementation seems to be for lagrangian and moment-based dof bases only. This code will not work if a developer implements a new Dof basis type. It seems that the mechanism of pushing the dof basis to the physical space should be in the abstract interface for dof bases so that it can be implemented also for new types. You are allocating temporary arrays. This part of the code is quite performance critical and allocating temporary arrays should be avoided. It looks that you are also evaluating the reference shape functions over and over e.g. /~https://github.com/Balaje/GSoC-2021/blob/f52c2a04a9367cac8c191ac287b167ad6a1dd9b2/Interpolation/interpolable.jl#L49 |
|
Hi @fverdugo We already discussed @Balaje @ericneiva and myself the point about how to generalise the implementation. We have observed, as you mention, that what we need to do is to implement @Balaje when you finish the part ☝️ we can start working on an efficient implementation. We can discuss that next week or via Slack. |
|
Hi @Balaje Here you have the efficient implementation of Please add this into the code, write some tests, and include it in this PR or open a separate one. Disclaimer. I haven't run the code. It can have typos. # Put this in ReferenceFEs/Dofs.jl
struct PushDofMap <: Map end
function evaluate!(cache,::PushDofMap,f::Dof,m::Field)
@abstractmethod
end
function evaluate!(cache,::PushDofMap,f::AbstractArray{<:Dof},m::Field)
@abstractmethod
end
# Local implementations
function replace_nodes(f::LagrangianDofBasis,x)
LagrangianDofBasis(x, f.dof_to_node, f.dof_to_comp, f.node_and_comp_to_dof)
end
function return_cache(::PushDofMap,f::LagrangianDofBasis,m::Field)
q = f.nodes
return_cache(m,q)
end
function evaluate!(cache,::PushDofMap,f::LagrangianDofBasis,m::Field)
q = f.nodes
x = evaluate!(cache,m,q)
replace_nodes(f,x)
end
function replace_nodes(f::MomentBasedDofBasis,x)
MomentBasedDofBasis(x, f.face_moments, f.face_nodes)
end
function return_cache(::PushDofMap,f::MomentBasedDofBasis,m::Field)
q = f.nodes
return_cache(m,q)
end
function evaluate!(cache,::PushDofMap,f::MomentBasedDofBasis,m::Field)
q = f.nodes
x = evaluate!(cache,m,q)
replace_nodes(f,x)
end
# For performance reasons we also want these global implementations
# in practice, only the global implementations will be used
function lazy_map(
::PushDofMap,
cell_f::AbstractArray{<:LagrangianDofBasis},
cell_m::AbstractArray{<:Field})
cell_q = lazy_map(f->f.nodes,cell_f)
cell_x = lazy_map(evaluate,cell_m,cell_q)
lazy_map(replace_nodes,cell_f,cell_x)
end
function lazy_map(
::PushDofMap,
cell_f::AbstractArray{<:MomentBasedDofBasis},
cell_m::AbstractArray{<:Field})
cell_q = lazy_map(f->f.nodes,cell_f)
cell_x = lazy_map(evaluate,cell_m,cell_q)
lazy_map(replace_nodes,cell_f,cell_x)
end
# In CellData/CellDofs.jl
# Replace this
function change_domain(a::CellDof,::ReferenceDomain,::PhysicalDomain)
@notimplemented
end
# by
function change_domain(a::CellDof,::ReferenceDomain,::PhysicalDomain)
trian = get_triangulation(a)
cell_m = get_cell_map(trian)
cell_f_ref = get_data(a)
cell_f_phys = lazy_map(PushDofMap(),cell_f_ref,cell_m)
CellDof(cell_f_phys,trian,DomainStyle(a))
end
|
|
After adding this you should be able to call |
|
Thanks, @fverdugo. I just want to point out that |
Once I get the basis s = change_domain(get_fe_dof_basis(new_fe_space), ReferenceDomain(), PhysicalDomain()) I do this inside basis = get_data(s)
cache = return_cache(testitem(bs), fh)
cell_vals = lazy_map(i-> evaluate!(cache, basis[i], fh), 1:num_cells(trian_new_space))Is getting Old method: /~https://github.com/Balaje/GSoC-2021/blob/main/Interpolation/interpolable_2.jl When I check julia> include("interpolable.jl")
julia> @btime FESpaces._cell_vals(V₂, Interpolatable(fh));
164.250 μs (4369 allocations: 449.45 KiB)But when I check julia> include("interpolable_2.jl");
julia> @btime FESpaces._cell_vals(V₂, Interpolatable(fh));
153.834 μs (3913 allocations: 416.05 KiB)I can see slightly better memory allocation and execution time in the old method. Does finding cell_vals = lazy_map(i-> evaluate!(cache, basis[i], fh), 1:num_cells(trian_new_space))create any problem? |
|
Here you using |
|
another issue is that |
|
Hi everyone, I pushed the changes with @fverdugo's interface along with the new Gridap.jl/src/CellData/CellFields.jl Line 861 in 7c3b3dc As @fverdugo mentioned, we could pass optional parameters like I now define the Gridap.jl/src/FESpaces/SingleFieldFESpaces.jl Line 310 in 7c3b3dc I checked the run-time and memory for the new implementation and I also compared the time it would take for # V₁_model = CartesianDiscreteModel((0,1,0,1),(2,2))
julia> fh = interpolate_everywhere(f, V₁);
julia> @btime FESpaces._cell_vals(V₁, fh); # Same Triangulation
5.597 μs (166 allocations: 17.59 KiB)
julia> ifh = Interpolable(fh);
# V₂_model = CartesianDiscreteModel((0,1,0,1),(4,4))
julia> @btime FESpaces._cell_vals(V₂, ifh);
16.041 μs (286 allocations: 31.48 KiB)
julia> gh = interpolate_everywhere(ifh, V₂); # Different Triangulation using the new interface
julia> @btime FESpaces._cell_vals(V₂, gh); # Same Triangulation
8.167 μs (216 allocations: 22.91 KiB)I also added the tests for |
There was a problem hiding this comment.
Hi @Balaje !
yes this is definitively better. I see a potential pitfall though. The array you return in _cell_vals cannot be consumed safely in a threaded loop. Since all threads would consume the same cache object stored in the Interpolatable
src/CellData/CellFields.jl
Outdated
| return_cache(f::Interpolable, x::Point) = f.cache | ||
| return_cache(f::Interpolable, xs::AbstractVector{<:Point}) = f.cache | ||
| evaluate!(cache, f::Interpolable, x::Point) = evaluate!(cache, f.uh, x) | ||
| evaluate!(cache, f::Interpolable, xs::AbstractVector{<:Point}) = evaluate!(cache, f.uh, xs) |
There was a problem hiding this comment.
Try this:
cache = return_cache(f,x)
evaluate!(cache,x)
# and
cache = return_cache(f,[x])
evaluate!(cache,[x])with f::Interpolatable and x::Point. I think one of the two, will fail.
There was a problem hiding this comment.
Hi @fverdugo Did you mean
evaluate!(cache, f, x)If so, it seems to work.
julia> x = Point(rand(2));
julia> cache = return_cache(ifh, x);
julia> evaluate!(cache, ifh, x)
1.0107218313827602
julia> cache = return_cache(ifh, [x]);
julia> evaluate!(cache, ifh, [x])
1-element Vector{Float64}:
1.0107218313827602
julia> typeof(ifh)
Interpolable{SingleFieldFEFunction{GenericCellField{ReferenceDomain}}}There was a problem hiding this comment.
ok! forget about this. this will change anyway by my comments below.
src/FESpaces/FESpaces.jl
Outdated
| @@ -181,6 +181,8 @@ export FESpaceWithLinearConstraints | |||
|
|
|||
| export FiniteElements | |||
|
|
|||
| export interpolate_everywhere_non_compatible_trian | |||
|
|
|||
src/FESpaces/SingleFieldFESpaces.jl
Outdated
| trian = get_triangulation(V) | ||
| cache = return_cache(f, Point(rand(2))) | ||
| b = change_domain(fe_basis, ReferenceDomain(), PhysicalDomain()) | ||
| cf = CellField(x->evaluate!(cache, f, x), trian, ReferenceDomain()) |
There was a problem hiding this comment.
This should be PhysicalDomain right?
even though I think has no effect since you are directly extracting the data get_data(cf)
|
I think I have it. Once the functionality to change the domain style of the dof basis is available, the implementation of struct KDTreeSearch end
struct Interpolable{M,A} <: Function
uh::A
tol::Float64
searchmethod::M
function Interpolable(uh; tol=1e-6, searchmethod=KDTreeSearch())
new{typeof(searchmethod),typeof(uh)}(uh, tol,searchmethod)
end
end
# Instead of delegating to uh (which takes default configuration for search method),
# one could implement these functions by using the extra parameters.
# One could even dispatch on search method (the type parameter M)
# In this case it is better that the search method is an instance
# of some struct instead a symbol.
# The implementation should look like the one for uh
return_cache(a::Interpolable,x::Point) = return_cache(a.uh,x)
return_cache(a::Interpolable,x::AbstractVector{<:Point}) = return_cache(a.uh,x)
evaluate!(cache,a::Interpolable,x::Point) = evaluate!(cache,a.uh,x)
evaluate!(cache,a::Interpolable,x::AbstractVector{<:Point}) = evaluate!(cache,a.uh,x)
function _cell_vals(V::SingleFieldFESpace, f::Interpolable)
σ = get_fe_dof_basis(V)
trian = get_triangulation(σ)
cell_f = Fill(GenericField(f),num_cells(trian))
cf = CellField(cell_f,trian,PhysicalDomain())
σ(cf)
end
This also should be thread safe. @Balaje can you try if this works ? |
|
And perhaps even further simplify to function _cell_vals(V::SingleFieldFESpace, f::Interpolable)
σ = get_fe_dof_basis(V)
trian = get_triangulation(σ)
cf = CellField(f,trian)
σ(cf)
end |
|
Even more simple: Perhaps you don't need to overload The key is just to define |
@fverdugo julia> fh
SingleFieldFEFunction():
num_cells: 4
DomainStyle: ReferenceDomain()
Triangulation: CartesianGrid()
Triangulation id: 13773240885803307439
julia> ifh = Interpolable(fh);
julia> @btime FESpaces._cell_vals(V₂, ifh);
245.083 μs (4582 allocations: 447.56 KiB)
julia> @btime FESpaces._cell_vals(V₁, fh);
5.757 μs (166 allocations: 17.59 KiB)
julia> @btime FESpaces._cell_vals(V₂, gh);
8.486 μs (216 allocations: 22.91 KiB) |
more than the absolute value, we want to see how this scales when refining the mesh |
|
Hi @fverdugo
I plot the memory allocated against the number of partition on each axis for a 2D mesh. Code Here. I observe that the allocations increase as I decrease the mesh size (the source mesh, where we build
With the new implementation, aren't we building the cache repeatedly while evaluating the
Also, can you explain this bit in detail as to how the old implementation is thread unsafe? Thanks |
yes, the cache stored inside the Interpolatable acts like a global variable. You can have race conditions when threads consume it. |
|
|
|
Hi @Balaje I believe that the elegant way of implementing this PR is just by defining the following code (in particular no overload of _cell_vals and it is even not necessary to implement the change_domain for DellDof) struct KDTreeSearch end
struct Interpolable{M,A} <: Function
uh::A
tol::Float64
searchmethod::M
function Interpolable(uh; tol=1e-6, searchmethod=KDTreeSearch())
new{typeof(searchmethod),typeof(uh)}(uh, tol,searchmethod)
end
end
Arrays.return_cache(a::Interpolable,x::Point) = return_cache(a.uh,x)
Arrays.return_cache(a::Interpolable,x::AbstractVector{<:Point}) = return_cache(a.uh,x)
Arrays.evaluate!(cache,a::Interpolable,x::Point) = evaluate!(cache,a.uh,x)
Arrays.evaluate!(cache,a::Interpolable,x::AbstractVector{<:Point}) = evaluate!(cache,a.uh,x)If the performance of this is not good, It means that some other parts of the code are not working properly. E.g., I have been playing a bit and I have found that evaluating a FEFunction on a point or at a vector of points allocats memory even with a cache: using Gridap
using Gridap.Arrays
using BenchmarkTools
k = 1
f(x) = x[1] + x[2]
reffe = ReferenceFE(lagrangian,Float64,k)
n = 10
model = CartesianDiscreteModel((0,1,0,1),(n,n))
V = FESpace(model,reffe)
vh = interpolate(f,V)
x = Point(.5,.45)
cache = return_cache(vh,x)
@btime evaluate!($cache,$vh,$x)
xs = [Point(.5,.45)]
cache = return_cache(vh,xs)
@btime evaluate!($cache,$vh,$xs)
This should be fixed otherwise the interpolation will not be efficient. Perhaps you can clean up the PR (i.e., only include the first code snippet in this comment and some tests) and in a second one we can fix the performance hot spots. |
@fverdugo Done.
I have been thinking about this. We could preallocate arrays, e.g., point_to_cell = map(x_to_cell,point_to_x)
cell_to_points,point_to_lpoint = make_inverse_table(point_to_cell,ncells)and store it in the cache. Is this a good place to start? Gridap.jl/src/CellData/CellFields.jl Line 364 in 8522d93 |
Hi @fverdugo @ericneiva @santiagobadia @oriolcg
Thanks for all your help on this. This is the PR for adding the features and tests for interpolation. It contains the following:
nlsolveto Newton Raphson and Backtracking to handle distorted meshes a little better.evaluate!forLagrangianDofBasis.evaluate!forMomentBasedDofBasis._cell_valsto generate cell-wise dof values.CellFieldto generate theGenericCellFieldobject for non-compatible triangulations.interpolate_everywhere_non_compatible_trianto generate theFEFunction. Also implemented for::MultiFieldFESpace.The only thing @ericneiva mentioned was whether to have
evaluate!(cache, b::LagrangianDofBasis, field, points)separately from the original implementation ofevaluate!(cache, b::LagrangianDofBasis, field)or havepointsas an optional argument.The PR is not urgent @fverdugo @ericneiva. Please have a look at it whenever you are free. Thanks 🙂