Random “Point source is outside the mesh” error

[497e0bdf29873]

The following MWE randomly fails with a “Point source is outside the mesh” error in scifem.PointSource. Sometimes when you run it, it finishes without error, but more frequently, it fails.

import numpy as np
import scifem
from dolfinx import fem, mesh
from mpi4py import MPI

domain = mesh.create_rectangle(
    MPI.COMM_WORLD,
    [
        np.array([-2, -2]),
        np.array([2, 2]),
    ],
    [32, 32],
    mesh.CellType.triangle,
)

V = fem.functionspace(domain, ("Lagrange", 2))

μ = [
    (np.array([0.3, 0.5]), 1.0),
    (np.array([0.7, 0.8]), 1.0),
    (np.array([1.0, 0.3]), 1.0),
]

for point, alpha in μ:
    ps = scifem.PointSource(V, point, magnitude=alpha)

The exact error:

$ python fubar.py
Traceback (most recent call last):
  File "fubar.py", line 25, in <module>
    ps = scifem.PointSource(V, point, magnitude=alpha)
  File "/opt/homebrew/Caskroom/miniconda/base/envs/fenicsx-env/lib/python3.13/site-packages/scifem/point_source.py", line 61, in __init__
    self.recompute_sources()
    ~~~~~~~~~~~~~~~~~~~~~~^^
  File "/opt/homebrew/Caskroom/miniconda/base/envs/fenicsx-env/lib/python3.13/site-packages/scifem/point_source.py", line 96, in recompute_sources
    raise ValueError("Point source is outside the mesh.")
ValueError: Point source is outside the mesh.

This is with scifem 0.13.1 and dolfinx 0.9.0. I tried updating everything to the latest versions available in conda, but it did not fix this error, and I ran into other problems with nanobind+dolfinx.

[jorgensd]

I was able to reproduce the error with a fresh conda installation (of 0.9 and 0.13.1).
However, with a local installation of scifem I can't reproduce, making it hard to debug.
(I also couldn't reproduce it in docker).

I also tried to upgrade the conda installation to 0.14.0, but that doesn't work either.

In your MWE, you are placing the point source at vertices of your mesh. Is that deliberate? If you know that the source is at a vertex, one can now use the vertex-measure dP in DOLFINx (since v0.10.0, initially added by @schnellerhase and extended by me) to specify integration over such entities.
Following is an MWE:

"""Create a point source aligning with a vertex in a distributed mesh.

Author: Jørgen S. Dokken
SPDX-License-Identifier: MIT
"""

from mpi4py import MPI
import dolfinx.fem.petsc
import ufl
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 5, 5)
V = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))

v = ufl.TestFunction(V)


# Find the relevant vertex in the distributed mesh and only add it on process that owns it
local_vertex = dolfinx.mesh.locate_entities(
    mesh, 0, lambda x: np.isclose(x[0], 0.2) & np.isclose(x[1], 0.4)
)
assert len(local_vertex) < 2
num_vertices_local = mesh.topology.index_map(0).size_local
if len(local_vertex) != 0 and local_vertex[0] < num_vertices_local:
    lv = local_vertex[0]
    mesh.topology.create_connectivity(0, mesh.topology.dim)
    v_to_c = mesh.topology.connectivity(0, mesh.topology.dim)
    c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)

    cell = v_to_c.links(lv)[0]
    local_vertex_idx = np.flatnonzero(c_to_v.links(cell) == lv)[0]
    itg_data = np.array([cell, local_vertex_idx], dtype=np.int32)
else:
    itg_data = np.empty(0, dtype=np.int32)


# Create integration measure for said pointsource
point_id = 33
itg_entity = [(point_id, itg_data)]
dP = ufl.Measure("dP", domain=mesh, subdomain_data=itg_entity, subdomain_id=point_id)

# Assemble RHS vector
magnitude = dolfinx.fem.Constant(mesh, 5.0)
L = magnitude * v * dP
b = dolfinx.fem.assemble_vector(dolfinx.fem.form(L))
b.scatter_reverse(dolfinx.la.InsertMode.add)
print(b.array)

if len(local_vertex) != 0 and local_vertex[0] < num_vertices_local:
    local_dof_idx = np.flatnonzero(b.array)
    dof_coords = V.tabulate_dof_coordinates()
    dof_coord = dof_coords[local_dof_idx[0]]
    print("Dof coordinate:", dof_coord)

[497e0bdf29873]

Appears I also cannot reproduce the problem on my Armbian box that has native packages available (but typical to centralised Linux package management, obsolete: scifem 0.10.0). Unfortunately, on macOS, Conda seems to be the only supported low-labour way to install Dolfinx. I wonder, though, if this problem has more to do with macOS or with Conda?

On the Armbian box I get many 'Authorization required, but no authorization protocol specified' complaints, which I don't get on macOS/Conda, so maybe this has something to do with how MPI is set up on each system?

The point sources matching vertices is not intentional. In fact, in the full code, they will be generated by an optimisation method.

[jorgensd]

Could you try v0.15.0: https://anaconda.org/channels/conda-forge/packages/scifem/overview and increase the new tol input to point source to say 1e-10?

[497e0bdf29873]

Unfortunately, doesn't help. Not even with tol=1e-1.

But I noticed something else odd. Inserting into recompute_sources the line

print(self._input_points, collision_data.dest_points, collision_data.dest_cells, collision_data.src_owner)

and running the script normally, I just get

[0.31 0.51] [] [] [-1]

If I mpirun -n 2, I get

[0.31 0.51] [] [] [0]
[0.31 0.51] [[3.10000000e-001 5.10000000e-001 2.83515055e-314]] [896] [-1]

If I mpirun -n 4, I get

[0.31 0.51] [] [] [-1]
[0.31 0.51] [] [] [-1]
[0.31 0.51] [] [] [-1]
[0.31 0.51] [[3.10000000e-001 5.10000000e-001 2.66107005e-314]] [246] [0]

So it seems that with MPI,

1. even workers that find the cell, fail, by returning src_owner=[-1], but also
2. Some MPI propagation still needs to be done? dolfinx.cpp.geometry.determine_point_ownership does not propagate?

(Not that I need MPI, I actually want to disable it as far as possible, and force a single worker, as how Dolfinx uses it, is not compatible with programs with a single controlling worker.)

[jorgensd]

So one thing to note is the documentation of PointSource:

scifem/src/scifem/point_source.py

Lines 41 to 44 in 4685c0e

	Note:
	Points should only be defined on one process. If they are sent in
	from multiple processes, multiple point sources will be created.

Points should only be defined on one process. If they are sent in from multiple processes, multiple point sources will be created.

It is of course strange that it finds the point (which are the same on all processes, with one search but not the others).

[jorgensd]

Could you try padding your points to be 3D?, i.e.

μ = [
    (np.array([0.3, 0.5, 0]), 1.0),
    (np.array([0.7, 0.8, 0]), 1.0),
    (np.array([1.0, 0.3, 0]), 1.0),
]

That resolves the issue for me.

[497e0bdf29873]

Yes! That solves it. It's the fixed 3-dimensionality of the C++-side, and apparently no safety/consistency checks there. Thanks!

[jorgensd]

I’ll add some simple checks to scifem to check that the points are 3D.