From 72012d729006b7d5144241e813594406d1c93f24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 9 Oct 2026 17:14:55 +0200 Subject: [PATCH 1/4] Fix NaN in parallel point evaluation with repeated knots At the right bound of the local subdomain, `TensorFemSpace.eval_field` and `eval_field_gradient` moved back by one knot from the span given by `find_span`, which selects an empty knot span (and yields NaN) if the bound is a knot of multiplicity larger than one. Add the option `side='left'` to `psydac.core.bsplines.find_span`, which selects the non-empty knot span on the left of an interior knot of any multiplicity. Use it in the new private method `TensorFemSpace._find_local_span`, which replaces the code duplicated in `eval_field` and `eval_field_gradient`. Add unit tests for `find_span` with `side='left'`, and an MPI test which evaluates a linear field and its gradient on the closed local subdomains with simple, repeated and mixed knot multiplicities. Fixes #612. --- psydac/core/bsplines.py | 19 +++++- psydac/core/tests/test_bsplines.py | 29 +++++++++ psydac/fem/tensor.py | 98 ++++++++++++++++-------------- psydac/fem/tests/test_tensor.py | 42 +++++++++++++ 4 files changed, 140 insertions(+), 48 deletions(-) diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 781a81626..64c0fb957 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -57,7 +57,7 @@ #============================================================================== -def find_span(knots, degree, x): +def find_span(knots, degree, x, side='right'): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -76,14 +76,29 @@ def find_span(knots, degree, x): x : float Location of interest. + side : {'right', 'left'}, default='right' + Which non-empty knot span to select if `x` is an interior knot, of any + multiplicity: the one on the right of `x`, such that + ``knots[span] <= x < knots[span+1]``, or the one on the left of `x`, + such that ``knots[span] < x <= knots[span+1]``. At the boundaries of the + domain the first or last knot span is always selected. + Returns ------- span : int Knot span index. """ + assert side in ('right', 'left') x = float(x) knots = np.ascontiguousarray(knots, dtype=float) - return find_span_p(knots, degree, x) + span = find_span_p(knots, degree, x) + + # Skip back over all the copies of x, if it is an interior knot + if side == 'left': + while span > degree and knots[span] == x: + span -= 1 + + return span #============================================================================== def find_spans(knots, degree, x, out=None): diff --git a/psydac/core/tests/test_bsplines.py b/psydac/core/tests/test_bsplines.py index d3cf20162..19200067e 100644 --- a/psydac/core/tests/test_bsplines.py +++ b/psydac/core/tests/test_bsplines.py @@ -39,6 +39,35 @@ def test_find_span( lims, nc, p, eps=1e-12 ): assert find_span( knots, p, x=xi ) == p + min( i, nc-1 ) assert find_span( knots, p, x=xi+eps ) == p + min( i, nc-1 ) +#============================================================================== +@pytest.mark.parametrize('p', (1, 2, 3)) +@pytest.mark.parametrize('multiplicity', ('simple', 'repeated', 'mixed')) +def test_find_span_side(p: int, multiplicity: str) -> None: + + grid = np.linspace(0.0, 1.0, 9) + interior = grid[1:-1] + m = {'simple' : [1] * len(interior), + 'repeated': [p] * len(interior), + 'mixed' : [1 + i % p for i in range(len(interior))]}[multiplicity] + knots = np.r_[[grid[0]] * (p + 1), np.repeat(interior, m), [grid[-1]] * (p + 1)] + first_span, last_span = p, len(knots) - p - 2 + + # On the breakpoints, the selected span is never empty + for x in grid: + span_right = find_span(knots, p, x, side='right') + span_left = find_span(knots, p, x, side='left') + if x == grid[0]: + assert span_right == span_left == first_span + elif x == grid[-1]: + assert span_right == span_left == last_span + else: + assert knots[span_right] <= x < knots[span_right + 1] + assert knots[span_left] < x <= knots[span_left + 1] + + # Away from the breakpoints, the side does not matter + for x in (grid[:-1] + grid[1:]) / 2: + assert find_span(knots, p, x, side='left') == find_span(knots, p, x) + #============================================================================== @pytest.mark.parametrize( 'lims', ([0,1], [-2,3]) ) @pytest.mark.parametrize( 'nc', (10, 18, 33) ) diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 527e6f50c..31c292bdc 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -210,6 +210,49 @@ def is_vector_valued(self): #-------------------------------------------------------------------------- # Abstract interface: evaluation methods #-------------------------------------------------------------------------- + @staticmethod + def _find_local_span(x: float, xlim: tuple, space: SplineSpace) -> tuple: + """ + Find the knot span of a point, as seen from the local interval of the + process along one direction. + + Parameters + ---------- + x : float + The coordinate of the point along one direction. + + xlim : tuple + The bounds of the local interval along the same direction. + + space : SplineSpace + The 1D spline space along the same direction. + + Returns + ------- + x_local : float + The coordinate `x`, moved onto the closest bound of `xlim` if it lies + outside of it by at most one machine epsilon. + + span : int + The index of the knot span which contains `x_local` within `xlim`. + """ + scale = max(abs(xlim[0]), abs(xlim[1]), 1.0) + tol = np.spacing(scale) + + if x < xlim[0] and xlim[0] - x <= tol: + x_local = xlim[0] + elif x > xlim[1] and x - xlim[1] <= tol: + x_local = xlim[1] + else: + x_local = x + + # On the right bound of the local interval, the knot span on the right + # of x belongs to the next process: use the one on the left of x + side = 'left' if x_local == xlim[1] else 'right' + span = find_span(space.knots, space.degree, x_local, side=side) + + return x_local, span + def eval_field( self, field, *eta, weights=None): assert isinstance( field, FemField ) @@ -226,36 +269,17 @@ def eval_field( self, field, *eta, weights=None): field.coeffs.update_ghost_regions() for (x, xlim, space) in zip( eta, self.eta_lims, self.spaces ): - scale = max(abs(xlim[0]), abs(xlim[1]), 1.0) - tol = np.spacing(scale) - - if x < xlim[0] and xlim[0] - x <= tol: - x = xlim[0] - - elif x > xlim[1] and x - xlim[1] <= tol: - x = xlim[1] - - knots = space.knots - degree = space.degree - span = find_span( knots, degree, x ) - - #-------------------------------------------------# - # Fix span for boundaries between subdomains # - #-------------------------------------------------# - # TODO: Use local knot sequence instead of global # - # one to get correct span in all situations # - #-------------------------------------------------# - if x == xlim[1] and x != knots[-1-degree]: - span -= 1 - #-------------------------------------------------# - basis = basis_funs( knots, degree, x, span) + x_local, span = self._find_local_span(x, xlim, space) + knots = space.knots + degree = space.degree + basis = basis_funs( knots, degree, x_local, span) # If needed, rescale B-splines to get M-splines if space.basis == 'M': basis *= space.scaling_array[span-degree : span+1] # Determine local span - wrap_x = space.periodic and x > xlim[1] + wrap_x = space.periodic and x_local > xlim[1] loc_span = span - space.nbasis if wrap_x else span bases.append( basis ) @@ -635,29 +659,11 @@ def eval_field_gradient(self, field, *eta, weights=None): index = [] for (x, xlim, space) in zip( eta, self.eta_lims, self.spaces ): - scale = max(abs(xlim[0]), abs(xlim[1]), 1.0) - tol = np.spacing(scale) - - if x < xlim[0] and xlim[0] - x <= tol: - x = xlim[0] - - elif x > xlim[1] and x - xlim[1] <= tol: - x = xlim[1] - + x_local, span = self._find_local_span(x, xlim, space) knots = space.knots degree = space.degree - span = find_span( knots, degree, x ) - #-------------------------------------------------# - # Fix span for boundaries between subdomains # - #-------------------------------------------------# - # TODO: Use local knot sequence instead of global # - # one to get correct span in all situations # - #-------------------------------------------------# - if x == xlim[1] and x != knots[-1-degree]: - span -= 1 - #-------------------------------------------------# - basis_0 = basis_funs(knots, degree, x, span) - basis_1 = basis_funs_1st_der(knots, degree, x, span) + basis_0 = basis_funs(knots, degree, x_local, span) + basis_1 = basis_funs_1st_der(knots, degree, x_local, span) # If needed, rescale B-splines to get M-splines if space.basis == 'M': @@ -666,7 +672,7 @@ def eval_field_gradient(self, field, *eta, weights=None): basis_1 *= scaling # Determine local span - wrap_x = space.periodic and x > xlim[1] + wrap_x = space.periodic and x_local > xlim[1] loc_span = span - space.nbasis if wrap_x else span bases_0.append( basis_0 ) diff --git a/psydac/fem/tests/test_tensor.py b/psydac/fem/tests/test_tensor.py index bf3610ef8..d6e5e00a6 100644 --- a/psydac/fem/tests/test_tensor.py +++ b/psydac/fem/tests/test_tensor.py @@ -5,6 +5,7 @@ #---------------------------------------------------------------------------# import os import contextlib +from itertools import product from pathlib import Path import pytest @@ -20,6 +21,10 @@ from psydac.mapping.discrete_gallery import discrete_mapping from psydac.api.discretization import discretize +from psydac.ddm.cart import DomainDecomposition +from psydac.fem.basic import FemField +from psydac.fem.splines import SplineSpace +from psydac.fem.tensor import TensorFemSpace #============================================================================== @@ -223,3 +228,40 @@ def test_plot_2d_decomposition(kind, root): test_plot_2d_decomposition('spline', 'first') test_plot_2d_decomposition('analytical', 'last') plt.show() + +#============================================================================== +# The right bound of each local subdomain (except the last) is an interior knot, +# whose multiplicity may be larger than one +@pytest.mark.mpi +@pytest.mark.parametrize('multiplicity', ['simple', 'repeated', 'mixed']) +def test_TensorFemSpace_eval_field_parallel(multiplicity: str) -> None: + + p = 3 + grid = np.linspace(0.0, 1.0, 9) + interior = grid[1:-1] + # Mixed: repeated knots at 0.25, 0.5 and 0.75, simple knots in between + m = {'simple' : [1] * len(interior), + 'repeated': [p] * len(interior), + 'mixed' : [1, 2, 1, 3, 1, 2, 1]}[multiplicity] + knots = np.r_[[grid[0]] * (p + 1), np.repeat(interior, m), [grid[-1]] * (p + 1)] + + spaces = [SplineSpace(degree=p, knots=knots) for _ in range(2)] + comm = MPI.COMM_WORLD + domain_decomposition = DomainDecomposition([8, 8], [False, False], comm=comm) + V = TensorFemSpace(domain_decomposition, *spaces) + + # Field f(x, y) = x + 2y, whose coefficients are given by the Greville abscissae + g1, g2 = [W.greville for W in V.spaces] + s1, s2 = V.coeff_space.starts + e1, e2 = V.coeff_space.ends + f = FemField(V) + f.coeffs[s1:e1+1, s2:e2+1] = g1[s1:e1+1, None] + 2 * g2[None, s2:e2+1] + f.coeffs.update_ghost_regions() + + # Evaluate on the closed local subdomain, including its bounds + starts, ends = V.local_domain + grids = [np.linspace(W.breaks[s], W.breaks[e + 1], 5) + for W, s, e in zip(V.spaces, starts, ends)] + for x, y in product(*grids): + assert np.isclose(f(x, y), x + 2 * y, rtol=0, atol=1e-14) + assert np.allclose(f.gradient(x, y), [1.0, 2.0], rtol=0, atol=1e-13) From d780a181a2b5fee36f0ef0953180c00c138ef8fd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 9 Oct 2026 17:53:07 +0200 Subject: [PATCH 2/4] Make eval_field_gradient behave like eval_field Unlike `TensorFemSpace.eval_field`, `eval_field_gradient` did not update the ghost regions of the field coefficients when they were not in sync, and did not check the space of the weights. Add both, and check the weights with `is not None` in both methods. Add helper functions to the tests in `psydac.fem.tests.test_tensor`, and add an MPI test which evaluates the field or its gradient (separately) on the closed local subdomains, starting from coefficients whose ghost regions are not in sync. --- psydac/fem/tensor.py | 12 ++++-- psydac/fem/tests/test_tensor.py | 66 +++++++++++++++++++++++++-------- 2 files changed, 59 insertions(+), 19 deletions(-) diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 31c292bdc..fe6c1d148 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -258,7 +258,7 @@ def eval_field( self, field, *eta, weights=None): assert isinstance( field, FemField ) assert field.space is self assert len( eta ) == self.ldim - if weights: + if weights is not None: assert weights.space == field.coeffs.space bases = [] @@ -287,7 +287,7 @@ def eval_field( self, field, *eta, weights=None): # Get contiguous copy of the spline coefficients required for evaluation index = tuple( index ) coeffs = field.coeffs[index].copy() - if weights: + if weights is not None: coeffs *= weights[index] # Evaluation of multi-dimensional spline @@ -653,11 +653,17 @@ def eval_field_gradient(self, field, *eta, weights=None): assert isinstance(field, FemField) assert field.space is self assert len(eta) == self.ldim + if weights is not None: + assert weights.space == field.coeffs.space bases_0 = [] bases_1 = [] index = [] + # Necessary if vector coeffs is distributed across processes + if not field.coeffs.ghost_regions_in_sync: + field.coeffs.update_ghost_regions() + for (x, xlim, space) in zip( eta, self.eta_lims, self.spaces ): x_local, span = self._find_local_span(x, xlim, space) knots = space.knots @@ -682,7 +688,7 @@ def eval_field_gradient(self, field, *eta, weights=None): # Get contiguous copy of the spline coefficients required for evaluation index = tuple( index ) coeffs = field.coeffs[index].copy() - if weights: + if weights is not None: coeffs *= weights[index] # Evaluate each component of the gradient using algorithm described in "Option 1" above diff --git a/psydac/fem/tests/test_tensor.py b/psydac/fem/tests/test_tensor.py index d6e5e00a6..6ba6b7161 100644 --- a/psydac/fem/tests/test_tensor.py +++ b/psydac/fem/tests/test_tensor.py @@ -229,6 +229,33 @@ def test_plot_2d_decomposition(kind, root): test_plot_2d_decomposition('analytical', 'last') plt.show() +#============================================================================== +def make_linear_field(knots: np.ndarray, degree: int) -> FemField: + """ + Create the field f(x, y) = x + 2y on a 2D spline space with the same knots + along both directions, distributed over all processes. The ghost regions + of its coefficients are not updated. + """ + spaces = [SplineSpace(degree=degree, knots=knots) for _ in range(2)] + ncells = [len(W.breaks) - 1 for W in spaces] + comm = MPI.COMM_WORLD + V = TensorFemSpace(DomainDecomposition(ncells, [False, False], comm=comm), *spaces) + + # The coefficients of a linear field are given by the Greville abscissae + g1, g2 = [W.greville for W in V.spaces] + s1, s2 = V.coeff_space.starts + e1, e2 = V.coeff_space.ends + f = FemField(V) + f.coeffs[s1:e1+1, s2:e2+1] = g1[s1:e1+1, None] + 2 * g2[None, s2:e2+1] + return f + +def local_points(V: TensorFemSpace) -> list: + """Return a grid of points on the closed local subdomain, including its bounds.""" + starts, ends = V.local_domain + grids = [np.linspace(W.breaks[s], W.breaks[e + 1], 5) + for W, s, e in zip(V.spaces, starts, ends)] + return list(product(*grids)) + #============================================================================== # The right bound of each local subdomain (except the last) is an interior knot, # whose multiplicity may be larger than one @@ -245,23 +272,30 @@ def test_TensorFemSpace_eval_field_parallel(multiplicity: str) -> None: 'mixed' : [1, 2, 1, 3, 1, 2, 1]}[multiplicity] knots = np.r_[[grid[0]] * (p + 1), np.repeat(interior, m), [grid[-1]] * (p + 1)] - spaces = [SplineSpace(degree=p, knots=knots) for _ in range(2)] - comm = MPI.COMM_WORLD - domain_decomposition = DomainDecomposition([8, 8], [False, False], comm=comm) - V = TensorFemSpace(domain_decomposition, *spaces) - - # Field f(x, y) = x + 2y, whose coefficients are given by the Greville abscissae - g1, g2 = [W.greville for W in V.spaces] - s1, s2 = V.coeff_space.starts - e1, e2 = V.coeff_space.ends - f = FemField(V) - f.coeffs[s1:e1+1, s2:e2+1] = g1[s1:e1+1, None] + 2 * g2[None, s2:e2+1] + f = make_linear_field(knots, p) f.coeffs.update_ghost_regions() - # Evaluate on the closed local subdomain, including its bounds - starts, ends = V.local_domain - grids = [np.linspace(W.breaks[s], W.breaks[e + 1], 5) - for W, s, e in zip(V.spaces, starts, ends)] - for x, y in product(*grids): + for x, y in local_points(f.space): assert np.isclose(f(x, y), x + 2 * y, rtol=0, atol=1e-14) assert np.allclose(f.gradient(x, y), [1.0, 2.0], rtol=0, atol=1e-13) + +#============================================================================== +# Both evaluation methods must update the ghost regions if they are not in sync +@pytest.mark.mpi +@pytest.mark.parametrize('method', ['field', 'gradient']) +def test_TensorFemSpace_eval_field_ghost_regions(method: str) -> None: + + p = 3 + grid = np.linspace(0.0, 1.0, 9) + knots = np.r_[[grid[0]] * p, grid, [grid[-1]] * p] + + f = make_linear_field(knots, p) + assert not f.coeffs.ghost_regions_in_sync + + for x, y in local_points(f.space): + if method == 'field': + assert np.isclose(f(x, y), x + 2 * y, rtol=0, atol=1e-14) + else: + assert np.allclose(f.gradient(x, y), [1.0, 2.0], rtol=0, atol=1e-13) + + assert f.coeffs.ghost_regions_in_sync From c9b71cc8381514cec51accb19b4420346559690d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Sat, 10 Oct 2026 00:19:35 +0200 Subject: [PATCH 3/4] Mention fixes of field evaluation in CHANGELOG.md --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index d8fc460e8..81360debd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ All notable changes to this project will be documented in this file. ### Added +- #614 : Add option `side` to `psydac.core.bsplines.find_span` to select the knot span on the left of a knot - #576 : Add module `psydac.feec.polar.conga_projections` with broken FEEC polar projections in 2D - #576 : Add 2D Poisson and TE Maxwell examples on polar mapped domains in `psydac/feec/polar/examples` - #576 : Add module `psydac.utilities.parallel_utils` for parallel execution and gathering variable-length arrays @@ -15,6 +16,8 @@ All notable changes to this project will be documented in this file. ### Fixed +- #614 : Fix NaN in parallel evaluation of fields at subdomain boundaries with repeated knots, see #612 +- #614 : Update ghost regions in `TensorFemSpace.eval_field_gradient` if needed, as in `eval_field` - #576 : Fix bug in `TensorFemSpace.eval_field` caused by round-off at MPI subdomain boundaries - #576 : Require `sympde==0.19.3` which fixes a bug in the linearity checks - #579 : Require `h5py>=3.16` which installs correctly with `setuptools>=81.0` From 40c9ef037de1900bd05e951951a4caa15759b4fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 19 Oct 2026 20:56:16 +0200 Subject: [PATCH 4/4] Remove unused import of os in tensor.py --- psydac/fem/tensor.py | 1 - 1 file changed, 1 deletion(-) diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index fe6c1d148..957611b34 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -8,7 +8,6 @@ of compact support """ -import os import itertools from types import MappingProxyType