Fix set_exterior_facet_domains for arbitrary ghosting

[jpdean]

Currently, set_exterior_facet_domains makes assumptions about the ghosting of the mesh. Looks like it was missed when #2337 was merged. This PR removes these assumptions. Unfortunately, it relies on an std::lower_bound in a loop, which is not ideal. However, I plan to replace this in #2269

[chrisrichardson]

This fixes an issue, and is no worse than what went before, in terms of efficiency (replaces binary search with lower bound)...

[francesco-ballarin]

@jpdean is this supposed to affect serial runs too, or not?

I have a few tutorials in multiphenicsx that worked before this PR was merged, and failed tonight, in some cases returning nan on solve, in other cases returning a value of dolfinx.fem.assemble_scalar which is not the expected one.

The failing tutorials are

tutorial_3a_advection_diffusion_reaction_neumann_control.ipynb
tutorial_4a_poisson_dirichlet_control.ipynb
tutorial_4b_poisson_neumann_control_boundary_observation.ipynb
tutorial_8_navier_stokes_neumann_control.ipynb

in
https://github.com/multiphenics/multiphenicsx/tree/main/tutorials/06_optimal_control
while

tutorial_3b_advection_diffusion_reaction_neumann_control.ipynb
tutorial_7a_stokes_dirichlet_control.ipynb
tutorial_7b_stokes_neumann_control.ipynb

are tutorials that use ds or dS too, but are not failing.

The only difference I can spot between the two groups is that in the failing ones I only define a custom ds measure and use ufl.dx for volume integration, while in the ones that do not fail I define either a custom dx or a custom dS measure too. Can this really be the relevant difference? Is there a unit test covering both cases?

[chrisrichardson]

I think there may be a lack of suitable unit tests. The PR should have fixed a bug in correctly determining surface facets for ds integrals. Would it be possible for you to create a minimal failing example from one of your notebooks?

[jpdean]

Oh no! Thanks for reporting this. I'll do some tests to see if I can figure out what is going wrong. As Chris mentioned, a minimum failing example would be really helpful.

[chrisrichardson]

@francesco-ballarin - it seems that the problem is due to unsorted MeshTags. The indices need to be sorted (it is stated in the documentation: https://docs.fenicsproject.org/dolfinx/main/cpp/mesh.html#_CPPv4I0EN7dolfinx4mesh8MeshTagsE)
There is a test for this in Debug mode, but it is not checked in Release mode, so it has caused this problem. Maybe we should have the check permanently on, given that it can cause silent wrong answers...

[garth-wells]

I don't think we should have the check on by default - that's what we have DEBUG mode for. Release mode checks should be O(1).

[jorgensd]

I guess it should be added to the Python documentation, as most Python users will not dig into the C++ documentation.
https://docs.fenicsproject.org/dolfinx/v0.5.1/python/generated/dolfinx.mesh.html#dolfinx.mesh.meshtags

[francesco-ballarin]

Thanks @chrisrichardson for sorting (pun intended ;) ) this out! I've pushed a fix to my repo.

I've gone through git blame of both dolfinx and multiphenicsx and apparently the timeline was the following:

• on March 30, 2020, MeshTags were introduced in dolfinx, without any sorting requirement
• on April 9, 2020, I updated my tutorials to use MeshTags
• on April 11, 2020, the sorting requirement was added in dolfinx, and somehow I missed or overlooked that ;)
  This was clearly my fault, but I agree that it might be helpful to add the requirement to the python documentation too. Still, I am amazed that this mistake only came to bite me 2.5 years later! :)