providing marked_vertices for InterfaceCoupling decreases robustness

[richardschu]

Working on #433 I realized, that omitting the marked_vertices allowed finding all the points on the interface (using the test presented in #433). Providing the vertices is optional though, and the Docu from RPE notes the argument function marked_vertices as a
"Function that marks relevant vertices to make search of active cells around point more efficient."

get_marked_vertices_via_boundary_ids marks vertices on faces depending on some boundary ID, which are passed to InterfaceCoupling::setup().

I would argue that the "relevant vertices" are the ones on the whole cell -- depending on the implamentation details which I did not check (yet). Maybe the prarallel layout plays a role as well ... I will adapt the vertices in an alternative function and let you know.

[nfehn]

Interesting result! (From my current understanding, I would have expected the boundary vertices to be enough.) @kronbichler @peterrum maybe we could exchange ideas / discuss this topic in person, given that we had to adapt the functions finding (all) active cells around points in deal.II already in the past.

My opinion is that marked_vertices is for improved performance and may not deteriorate robustness. I think it should be ensured on the deal.II-side that robustness is maintained also when passing marked_vertices.

[richardschu]

I see the point that one might expect a "fallback"-option to a slower, but also more robust point search. However, I am also not sure if the input we provide is correct. In my little toy example the triangulation is a p::d::T and I tested with n=1 and n=4 cores, giving the following results:

i)
omitting the marked_vertices via
map_evaluator.emplace(quad_index, dealii::Utilities::MPI::RemotePointEvaluation<dim>(tolerance_, false, 0));
leads to a successful search (all points found for n=1 and n=4).
ii)
providing marked_vertices with all vertices on those faces with a certain boundary ID leads to 23 (n=4) or 24 (n=1) points not found
iii)
providing marked_vertices with all vertices of cells with any face on the boundary with the sought ID leads to 3 (n=4) or 4 (n=1) unmatched points ... so, already better, but I would expect this to give the same result as i).
iv)
providing a dummy marked_vertices marking all vertices in the tria (not improving performance but simply testing equivalence to the default version taken when no marked_vertices are provided) via
std::vector<bool> marked_vertices(dof_handler_src_.get_triangulation().n_vertices(), true); map_evaluator.emplace(quad_index, dealii::Utilities::MPI::RemotePointEvaluation<dim>( tolerance_, false, 0, [marked_vertices](){ return marked_vertices; }));
gives identical results as iii).

So the parallel layout seems to play a role (triagulation.n_vertices() returns the total number of vertices, also some maybe not used ones, see: )

So this could mean that:
a) we interpret the marked_vertices wrong (but option iv) is quite fail-proof isn't it?)
or
b) something is done differently in those 2 code paths inside deal.II which leads to the difference in results...

[peterrum]

To me it looks like that there are functions that use the vertices of the triangulation directly instead of the mapped ones, which would be returned by GridTools::Cache::get_used_vertices(). The function

https://github.com/dealii/dealii/blob/d47537fa246886e5308d5aecbf42358d73a00774/source/grid/grid_tools_dof_handlers.cc#L466-L468

calls

https://github.com/dealii/dealii/blob/d47537fa246886e5308d5aecbf42358d73a00774/source/grid/grid_tools_dof_handlers.cc#LL67C41-L67C49

There are other instances in the namespace GridTools that do not use the mapped points.

I guess why providing no marked_vertices and one filled with true gives different results is

https://github.com/dealii/dealii/blob/d47537fa246886e5308d5aecbf42358d73a00774/source/grid/grid_tools_dof_handlers.cc#L543-L547

[kronbichler]

I tried to figure out where exactly we are arriving and ended up in this function here: https://github.com/dealii/dealii/blob/d47537fa246886e5308d5aecbf42358d73a00774/source/grid/grid_tools.cc#L2927-L2945 where I could not find the function in grid_tools_dof_handler.cc that does not use the mapping indicated by @peterrum above. Can you confirm. I then had a look in dealii/source/grid/grid_tools.cc and it seems to contain several bugs regarding the use of the mapping. Let me collect my findings here, because I do not have the time to check a solution yet.

• This one computes a distance between a mapped point and the plain vertices (probably wrong):
  https://github.com/dealii/dealii/blob/d47537fa246886e5308d5aecbf42358d73a00774/source/grid/grid_tools.cc#L3078
• By the way, the use of std::map in https://github.com/dealii/dealii/blob/d47537fa246886e5308d5aecbf42358d73a00774/source/grid/grid_tools.cc#L6459-L6475 is bad because this is probably the most expensive part of the search operation, potentially of quadratic complexity.

I might post more if I see more later.

[peterrum]

Can you confirm.

I guess I pointed out the wrong find_active_cell_around_point() instance. There are just too many... Anyway, both files grid_tools.cc and grid_tools_dof_handlers.cc have the same problems.

[kronbichler]

I agree, there are too many functions with similar tasks and too many wrinkles in the code that jump back and forth. It's just that I do not see a simple solution, either. What I would like to do, in case we identify the problem, is to check if the overall complexities are ok or if there is something to be done. I fear that those functions do a lot of useless work sometimes, but it is hard to look at it without a specific hint.

[peterrum]

I think the easiest approach for testing would be to take MappingQCache, move the vertices according this, and reset MappingQCache (similar to flatting in the case of refinement). In this case, the vertices and the mapped points should be identical for linear mapping. If this is successful, the missing mapping is the issue.

@richardschu Does this make sense to you. Could you try this out?

[richardschu]

Thanks for the help everyone! I think we managed to resolve that problem. In all my tests, the discretization error in the map compared to the tolerance is the cause if RPE is not succesful in finding any points. This is expected behavior and with a PR I will push to deal.II soon, this problem is fixed.
Regarding the vertices to provide as a hint for RPE: as far as I can tell, the closest vertex to the point within the marked_vertices is used to get the start cell batch (the connected cells + neighbors). Removing the early loop break, the more general point search (same as if no marked_cells where provided) is started.
Bonus comment: I believe that providing the interface vertices as marked_vertices as done already in ExaDG is the right way to go (and not using all the vertices of cells with a face on the interface).

[richardschu]

So I suggest we close this issue if you agree?

[kronbichler]

Yes, this sounds good to me. I am relieved that things work as they should now (at least when the comments are fixed). Should we wait with closing until the relevant fixes are in deal.II?

[nfehn]

yes, let's wait until the problem is really resolved, i.e. dealii/master has been updated accordingly.

[nfehn]

this issue seems to be resolved now?