diff --git a/CHANGELOG.md b/CHANGELOG.md index 60f90a5ec..3f3fbe36f 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 - #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 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..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 @@ -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,36 +268,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 ) @@ -263,7 +286,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 @@ -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 diff --git a/psydac/fem/tests/test_tensor.py b/psydac/fem/tests/test_tensor.py index bf3610ef8..6ba6b7161 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,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