Fix NaN in parallel point evaluation with repeated knots

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
 -   #527 : Add constructor `Geometry.from_file` with option `domain` to reuse an existing SymPDE domain, after checking that it matches the file
 -   #527 : Add function `get_available_mappings` and option `mpi_dims_mask` to `discrete_mapping` in `psydac.mapping.discrete_gallery`
 -   #527 : Add property `mpi_dims_mask` to `DomainDecomposition`
@@ -25,6 +26,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`
 -   #624 : Terminate all MPI processes when one of them fails in `parallel_run_from_cli`, see #621
 -   #604 : Fix several bugs in iterative linear solvers (MINRES, LSMR, GMRES, CG, BiCGSTAB), see #603
 -   #527 : Require `sympde==0.20.0` which changes how multipatch interfaces are defined

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):

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) )

psydac/fem/tensor.py

@@ -8,7 +8,6 @@
 of compact support
 
 """
-import os
 import itertools
 from types import MappingProxyType
 
@@ -210,12 +209,55 @@ 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 )
         assert field.space is self
         assert len( eta ) == self.ldim
-        if weights:
+        if weights is not None:
             assert weights.space == field.coeffs.space
 
         bases = []
@@ -226,44 +268,25 @@ 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 )
             index.append( slice( loc_span-degree, loc_span+1 ) )
         # 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
@@ -629,35 +652,23 @@ 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   = []
 
-        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]
+        # 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
             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 +677,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 )
@@ -676,7 +687,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

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,74 @@ def test_plot_2d_decomposition(kind, root):
     test_plot_2d_decomposition('spline', 'first')
     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
+@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)]
+
+    f = make_linear_field(knots, p)
+    f.coeffs.update_ghost_regions()
+
+    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