Efficient alternative for locate_cell to find all elements that share a point

[SouzaEM]

Hello, I have a similar issue than #3015. I need the cell that contains point $x_s$ and, if $x_s$ is on any internal edge of the mesh, I need both cells that share the edge. Here is the context of my problem: I want to compute the linear functional:

$\int_\Omega M_{ij}\frac{\partial v_i}{\partial x_j} dx$

where $v_i$ is a test function and $M_{ij}$ is a point moment tensor source of the form:

$M_{ij}(t, x) = \overline{M}_{ij} w(t) \delta(x - x_s)$

where $\overline{M}_{ij}$ is tensor, $w$ is a wavelet, $\delta$ is Dirac's delta, and $x_s$ is the source location. Then, we have:

$\int_\Omega M_{ij}\frac{\partial v_i}{\partial x_j} dx = \overline{M}_{ij} w(t) \frac{\partial v_i}{\partial x_j}(x_s)$

However, if $x_s$ is on the edge of an element, the right-hand side must be evaluated twice because $\frac{\partial v_i}{\partial x_j}$ is discontinuous on the edges (Mulder (2020)). Is there a efficient way of finding both cells given $x_s$. Currently, our software spyro uses locate_cell for other type of sources that do not involve the gradient of test functions.

[dham]

There are two issues here. One is the point evaluation functional. The other is your query about points that are on the edges of cells.

For the point evaluation functional, it sounds to me like Spyro is attempting to hand code something that Firedrake already explicitly supports. The right way to do this is to use function spaces defined on vertex only meshes. The exact case is described in https://doi.org/10.5194/gmd-17-5369-2024, especially listing 2. The same approach is used for an FWI demo in: https://www.firedrakeproject.org/demos/full_waveform_inversion.py.html and the VertexOnlyMesh is described in some more detail here: https://www.firedrakeproject.org/point-evaluation.html#primary-api-interpolation-onto-a-vertex-only-mesh

The second point is what to do on the boundaries of cells. The short version is that you don't actually ever know that you are on the boundary of a cell, because the roundoff errors mean that this is always uncertain. Even if you did know, the difference in even a discontinuous quantity from one side of the face to another is truncation error, so if it's at a level that matters to your answer this indicates that you're under-resolved. Averaging won't get you out of that.

Strictly speaking, point evaluation of a discontinuous function isn't actually well-defined anywhere so you're committing a variational crime by doing that point evaluation at all, but the magnitude of the crime is at most the truncation error (because you could have projected the quantity into a continuous space), so simply ignoring it and evaluating on whichever of the cells the algorithm finds is probably reasonable for most cases.

[SouzaEM]

I joined the spyro team a few months ago, but I believe they implemented the point evaluation some years ago before firedrake supported vertex meshes. Now that we updated the code to have good performance on newer firedrake versions, we can evaluate moving the implementation to vertex only meshes.

That is nice to know that the error of evaluating the derivative on only one element will be small. It will simplify the implementation.

Thank you for the detailed explanations and for taking your time to help me!

[Ig-dolci]

@SouzaEM Spyro recently merged an FWI solver (with automatic differentiation for adjoints) using VertexOnlyMesh.

[SouzaEM]

Daiane, I will check it out as soon as I merge the elastic wave solver in the main branch of spyro. Thank you for the suggestion!