From 458fa28ffe8b027f125b377fda11899ec163eee7 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 6 Sep 2023 16:36:37 +0200 Subject: [PATCH 01/37] add a test with higher multipliciy in the DeRham sequence --- psydac/api/tests/test_api_feec_2d.py | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 1da94f125..0ae0acc7b 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -185,7 +185,7 @@ def update_plot(fig, t, x, y, field_h, field_ex): def run_maxwell_2d_TE(*, use_spline_mapping, eps, ncells, degree, periodic, Cp, nsteps, tend, - splitting_order, plot_interval, diagnostics_interval, tol, verbose): + splitting_order, plot_interval, diagnostics_interval, tol, verbose, mult=1): import os @@ -290,7 +290,7 @@ class CollelaMapping2D(Mapping): #-------------------------------------------------------------------------- if use_spline_mapping: domain_h = discretize(domain, filename=filename, comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h) + derham_h = discretize(derham, domain_h, multiplicity = [mult,mult]) periodic_list = mapping.get_callable_mapping().space.periodic degree_list = mapping.get_callable_mapping().space.degree @@ -311,7 +311,7 @@ class CollelaMapping2D(Mapping): else: # Discrete physical domain and discrete DeRham sequence domain_h = discretize(domain, ncells=[ncells, ncells], periodic=[periodic, periodic], comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=[degree, degree]) + derham_h = discretize(derham, domain_h, degree=[degree, degree], multiplicity = [mult,mult]) # Discrete bilinear forms a1_h = discretize(a1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL) @@ -710,6 +710,25 @@ def test_maxwell_2d_periodic(): assert abs(namespace['error_Ey'] - ref['error_Ey']) / ref['error_Ey'] <= TOL assert abs(namespace['error_Bz'] - ref['error_Bz']) / ref['error_Bz'] <= TOL +def test_maxwell_2d_periodic_multiplicity(): + + namespace = run_maxwell_2d_TE( + use_spline_mapping = False, + eps = 0.5, + ncells = 12, + degree = 3, + periodic = True, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + verbose = False, + mult = 2 + ) + def test_maxwell_2d_dirichlet(): From 27027fe4235310291675e02c9288eba32f161b38 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 27 Sep 2023 20:25:48 +0200 Subject: [PATCH 02/37] problem seems to be fixed for non-periodic problems --- psydac/api/tests/test_api_feec_2d.py | 8 +++++--- psydac/core/bsplines.py | 4 ++-- psydac/feec/global_projectors.py | 2 +- psydac/fem/splines.py | 5 +++-- psydac/linalg/direct_solvers.py | 2 ++ 5 files changed, 13 insertions(+), 8 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 0ae0acc7b..6cb6ef5d3 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -715,9 +715,9 @@ def test_maxwell_2d_periodic_multiplicity(): namespace = run_maxwell_2d_TE( use_spline_mapping = False, eps = 0.5, - ncells = 12, - degree = 3, - periodic = True, + ncells = 10, + degree = 5, + periodic = False, Cp = 0.5, nsteps = 1, tend = None, @@ -847,6 +847,8 @@ def test_maxwell_2d_dirichlet_par(): #============================================================================== if __name__ == '__main__': + test_maxwell_2d_periodic_multiplicity() + exit() import argparse parser = argparse.ArgumentParser( diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 4bbb9a217..d82690b5b 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -354,7 +354,7 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None): return out #============================================================================== -def histopolation_matrix(knots, degree, periodic, normalization, xgrid, check_boundary=True, out=None): +def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multiplicity=1, check_boundary=True, out=None): """Computes the histopolation matrix. If called with normalization='M', this uses M-splines instead of B-splines. @@ -418,7 +418,7 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, check_bo knots = np.ascontiguousarray(knots, dtype=float) xgrid = np.ascontiguousarray(xgrid, dtype=float) - elevated_knots = elevate_knots(knots, degree, periodic) + elevated_knots = elevate_knots(knots, degree, periodic, multiplicity=multiplicity) normalization = normalization == "M" diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 773df5cf3..73d9ec4a5 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -163,7 +163,7 @@ def __init__(self, space, nquads = None): solverblocks += [KroneckerLinearSolver(tensorspaces[i].vector_space, solvercells)] - dataslice = tuple(slice(p, -p) for p in tensorspaces[i].vector_space.pads) + dataslice = tuple(slice(p*m, -p*m) for p,m in zip(tensorspaces[i].vector_space.pads,tensorspaces[i].vector_space.shifts)) dofs[i] = rhsblocks[i]._data[dataslice] # finish arguments and create a lambda diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index e18d09c94..f39afdc3d 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -132,7 +132,7 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul self._breaks = grid self._ncells = len(grid) - 1 self._greville = greville(knots, degree, periodic) - self._ext_greville = greville(elevate_knots(knots, degree, periodic), degree+1, periodic) + self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic) self._scaling_array = scaling_array self._parent_multiplicity = parent_multiplicity @@ -212,7 +212,8 @@ def init_histopolation( self, dtype=float): degree = self.degree, periodic = self.periodic, normalization = self.basis, - xgrid = self.ext_greville + xgrid = self.ext_greville, + multiplicity = self._multiplicity ) self.hmat= imat if self.periodic: diff --git a/psydac/linalg/direct_solvers.py b/psydac/linalg/direct_solvers.py index db6a49c2a..63aa295e6 100644 --- a/psydac/linalg/direct_solvers.py +++ b/psydac/linalg/direct_solvers.py @@ -149,6 +149,8 @@ def __init__( self, spmat ): assert isinstance( spmat, spmatrix ) self._space = np.ndarray + print("top") + print(spmat.tocsc()) self._splu = splu( spmat.tocsc() ) #-------------------------------------- From 4bdba960c407ee6e8f6d0651456dbf45f16d23b1 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Thu, 28 Sep 2023 13:25:52 +0200 Subject: [PATCH 03/37] working on periodic case but not working yet --- psydac/api/discretization.py | 1 - psydac/api/tests/test_api_feec_2d.py | 21 ++++++++++++++++++++- psydac/core/bsplines.py | 1 + psydac/core/bsplines_kernels.py | 9 ++++++--- psydac/fem/splines.py | 1 + 5 files changed, 28 insertions(+), 5 deletions(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index e097d4621..7c34f15bd 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -368,7 +368,6 @@ def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, # Create uniform grid grids = [np.linspace(xmin, xmax, num=ne + 1) for xmin, xmax, ne in zip(min_coords, max_coords, ncells)] - # Create 1D finite element spaces and precompute quadrature data spaces[i] = [SplineSpace( p, multiplicity=m, grid=grid , periodic=P) for p,m,grid,P in zip(degree_i, multiplicity_i,grids, periodic)] else: diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 6cb6ef5d3..7069ecd2d 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -710,7 +710,7 @@ def test_maxwell_2d_periodic(): assert abs(namespace['error_Ey'] - ref['error_Ey']) / ref['error_Ey'] <= TOL assert abs(namespace['error_Bz'] - ref['error_Bz']) / ref['error_Bz'] <= TOL -def test_maxwell_2d_periodic_multiplicity(): +def test_maxwell_2d_multiplicity(): namespace = run_maxwell_2d_TE( use_spline_mapping = False, @@ -728,6 +728,25 @@ def test_maxwell_2d_periodic_multiplicity(): verbose = False, mult = 2 ) + +def test_maxwell_2d_periodic_multiplicity(): + + namespace = run_maxwell_2d_TE( + use_spline_mapping = False, + eps = 0.5, + ncells = 4, + degree = 2, + periodic = True, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + verbose = False, + mult =2 + ) def test_maxwell_2d_dirichlet(): diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index d82690b5b..04149c61e 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -339,6 +339,7 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None): """ knots = np.ascontiguousarray(knots, dtype=float) xgrid = np.ascontiguousarray(xgrid, dtype=float) + if out is None: nb = len(knots) - degree - 1 if periodic: diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index dce0d059b..cf24e5948 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -811,7 +811,7 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] It should be of the appropriate shape and dtype. """ ncells = len(breaks) - 1 - for i in range(1, ncells): + for i in range(0, ncells + 1): out[degree + 1 + (i - 1) * multiplicity:degree + 1 + i * multiplicity] = breaks[i] out[degree] = breaks[0] @@ -820,8 +820,11 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] if periodic: period = breaks[-1]-breaks[0] - out[:degree] = breaks[ncells - degree:ncells] - period - out[len(out) - degree:] = breaks[1:degree + 1] + period + #out[:degree] = breaks[ncells - degree:ncells] - period + #out[len(out) - degree:] = breaks[1:degree + 1] + period + + out[:degree + 1-multiplicity] = out[(ncells-1) * multiplicity+1:degree + (ncells-2) * multiplicity + 1] - period + out[len(out) - degree + multiplicity-1:] = out[degree+multiplicity:2*degree + 1] + period else: out[0:degree + 1] = breaks[0] out[len(out) - degree - 1:] = breaks[-1] diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index f39afdc3d..1bbe49ba1 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -133,6 +133,7 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul self._ncells = len(grid) - 1 self._greville = greville(knots, degree, periodic) self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic) + self._scaling_array = scaling_array self._parent_multiplicity = parent_multiplicity From 0bcb11c224ed5a3ac962bc7d4971039cd738c93b Mon Sep 17 00:00:00 2001 From: vcarlier Date: Thu, 28 Sep 2023 13:31:05 +0200 Subject: [PATCH 04/37] remove print --- psydac/linalg/direct_solvers.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/psydac/linalg/direct_solvers.py b/psydac/linalg/direct_solvers.py index 63aa295e6..db6a49c2a 100644 --- a/psydac/linalg/direct_solvers.py +++ b/psydac/linalg/direct_solvers.py @@ -149,8 +149,6 @@ def __init__( self, spmat ): assert isinstance( spmat, spmatrix ) self._space = np.ndarray - print("top") - print(spmat.tocsc()) self._splu = splu( spmat.tocsc() ) #-------------------------------------- From 54ab751f99a15f2b20b150d2323b65e9079a8fc0 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Thu, 28 Sep 2023 13:33:53 +0200 Subject: [PATCH 05/37] remove added spaces or line --- psydac/api/discretization.py | 1 + psydac/core/bsplines.py | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index 7c34f15bd..e097d4621 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -368,6 +368,7 @@ def discretize_space(V, domain_h, *, degree=None, multiplicity=None, knots=None, # Create uniform grid grids = [np.linspace(xmin, xmax, num=ne + 1) for xmin, xmax, ne in zip(min_coords, max_coords, ncells)] + # Create 1D finite element spaces and precompute quadrature data spaces[i] = [SplineSpace( p, multiplicity=m, grid=grid , periodic=P) for p,m,grid,P in zip(degree_i, multiplicity_i,grids, periodic)] else: diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 04149c61e..d82690b5b 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -339,7 +339,6 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None): """ knots = np.ascontiguousarray(knots, dtype=float) xgrid = np.ascontiguousarray(xgrid, dtype=float) - if out is None: nb = len(knots) - degree - 1 if periodic: From c2d97a845d37777b6d8d6b6752be5e5b7de278dd Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 29 Sep 2023 17:10:09 +0200 Subject: [PATCH 06/37] changing the construction of the knots sequence to unable higher multiplicity nodes, few fixes to make it work --- psydac/api/tests/test_api_feec_2d.py | 9 +++---- psydac/core/bsplines.py | 19 ++++++++------ psydac/core/bsplines_kernels.py | 37 +++++++++++++++------------- psydac/feec/pull_push.py | 2 -- psydac/fem/splines.py | 9 +++++-- psydac/fem/tensor.py | 8 +++--- 6 files changed, 47 insertions(+), 37 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 7069ecd2d..96817ec7a 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -464,7 +464,7 @@ def discrete_energies(self, e, b): push_2d_hcurl(E.fields[0], E.fields[1], x1i, x2j, F) Bz_values[i, j] = push_2d_l2(B, x1i, x2j, F) - + # Electric field, x component fig2 = plot_field_and_error(r'E^x', x, y, Ex_values, Ex_ex(0, x, y), *gridlines) fig2.show() @@ -605,10 +605,9 @@ def discrete_energies(self, e, b): Ex_values[i, j], Ey_values[i, j] = \ push_2d_hcurl(E.fields[0], E.fields[1], x1i, x2j, F) - Bz_values[i, j] = push_2d_l2(B, x1i, x2j, F) - # ... + # ... # Error at final time error_Ex = abs(Ex_ex(t, x, y) - Ex_values).max() error_Ey = abs(Ey_ex(t, x, y) - Ey_values).max() @@ -734,8 +733,8 @@ def test_maxwell_2d_periodic_multiplicity(): namespace = run_maxwell_2d_TE( use_spline_mapping = False, eps = 0.5, - ncells = 4, - degree = 2, + ncells = 3, + degree = 3, periodic = True, Cp = 0.5, nsteps = 1, diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index d82690b5b..4803c8bfa 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, multiplicity = 1): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -83,7 +83,8 @@ def find_span(knots, degree, x): """ x = float(x) knots = np.ascontiguousarray(knots, dtype=float) - return find_span_p(knots, degree, x) + multiplicity = int(multiplicity) + return find_span_p(knots, degree, x, multiplicity = multiplicity) #============================================================================== def find_spans(knots, degree, x, out=None): @@ -301,7 +302,7 @@ def basis_funs_all_ders(knots, degree, x, span, n, normalization='B', out=None): return out #============================================================================== -def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None): +def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, multiplicity = 1): """Computes the collocation matrix If called with normalization='M', this uses M-splines instead of B-splines. @@ -349,8 +350,9 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None): assert out.shape == ((xgrid.shape[0], nb)) and out.dtype == np.dtype('float') bool_normalization = normalization == "M" + multiplicity = int(multiplicity) - collocation_matrix_p(knots, degree, periodic, bool_normalization, xgrid, out) + collocation_matrix_p(knots, degree, periodic, bool_normalization, xgrid, out, multiplicity=multiplicity) return out #============================================================================== @@ -433,8 +435,8 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli else: assert out.shape == (len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1) assert out.dtype == np.dtype('float') - - histopolation_matrix_p(knots, degree, periodic, normalization, xgrid, check_boundary, elevated_knots, out) + multiplicity = int(multiplicity) + histopolation_matrix_p(knots, degree, periodic, normalization, xgrid, check_boundary, elevated_knots, out, multiplicity = multiplicity) return out #============================================================================== @@ -607,7 +609,10 @@ def make_knots(breaks, degree, periodic, multiplicity=1, out=None): breaks = np.ascontiguousarray(breaks, dtype=float) if out is None: - out = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) + if periodic : + out = np.zeros(multiplicity * len(breaks[1:]) + 1 + 2 * degree) + else : + out = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) else: assert out.shape == (multiplicity * len(breaks[1:-1]) + 2 + 2 * degree,) \ and out.dtype == np.dtype('float') diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index cf24e5948..06489f18f 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -5,7 +5,7 @@ import numpy as np # ============================================================================= -def find_span_p(knots: 'float[:]', degree: int, x: float): +def find_span_p(knots: 'float[:]', degree: int, x: float, multiplicity : int = 1): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -36,9 +36,10 @@ def find_span_p(knots: 'float[:]', degree: int, x: float): """ # Knot index at left/right boundary low = degree - high = len(knots)-1-degree + high = len(knots)-multiplicity-degree # Check if point is exactly on left/right boundary, or outside domain + if x <= knots[low ]: return low if x >= knots[high]: return high-1 @@ -55,7 +56,7 @@ def find_span_p(knots: 'float[:]', degree: int, x: float): # ============================================================================= -def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]'): +def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', multiplicity : int = 1): """ Determine the knot span index at a set of locations x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -81,7 +82,7 @@ def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]'): n = x.shape[0] for i in range(n): - out[i] = find_span_p(knots, degree, x[i]) + out[i] = find_span_p(knots, degree, x[i], multiplicity=multiplicity) # ============================================================================= @@ -383,7 +384,7 @@ def basis_integrals_p(knots: 'float[:]', degree: int, out: 'float[:]'): # ============================================================================= def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normalization: bool, xgrid: 'float[:]', - out: 'float[:,:]'): + out: 'float[:,:]', multiplicity : int = 1): """ Compute the collocation matrix :math:`C_ij = B_j(x_i)`, which contains the values of each B-spline basis function :math:`B_j` at all locations :math:`x_i`. @@ -421,11 +422,10 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali basis = np.zeros((nx, degree + 1)) spans = np.zeros(nx, dtype=int) - find_spans_p(knots, degree, xgrid, spans) + find_spans_p(knots, degree, xgrid, spans, multiplicity = multiplicity) basis_funs_array_p(knots, degree, xgrid, spans, basis) # Fill in non-zero matrix values - # Rescaling of B-splines, to get M-splines if needed if not normalization: if periodic: @@ -460,7 +460,7 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali # ============================================================================= def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normalization: bool, xgrid: 'float[:]', - check_boundary: bool, elevated_knots: 'float[:]', out: 'float[:,:]'): + check_boundary: bool, elevated_knots: 'float[:]', out: 'float[:,:]', multiplicity : int = 1): """Computes the histopolation matrix. If called with normalization=True, this uses M-splines instead of B-splines. @@ -545,7 +545,8 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma False, False, xgrid_new[:actual_len], - colloc) + colloc, + multiplicity = multiplicity) m = colloc.shape[0] - 1 n = colloc.shape[1] - 1 @@ -811,21 +812,23 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] It should be of the appropriate shape and dtype. """ ncells = len(breaks) - 1 - for i in range(0, ncells + 1): - out[degree + 1 + (i - 1) * multiplicity:degree + 1 + i * multiplicity] = breaks[i] - - out[degree] = breaks[0] - out[len(out) - degree - 1] = breaks[-1] - + if periodic: + for i in range(1, ncells+1): + out[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] period = breaks[-1]-breaks[0] #out[:degree] = breaks[ncells - degree:ncells] - period #out[len(out) - degree:] = breaks[1:degree + 1] + period + out[:degree + 1] = out[ncells * multiplicity : ncells * multiplicity + degree + 1] - period + out[len(out) - degree :] = out[degree+1:2*degree + 1] + period + - out[:degree + 1-multiplicity] = out[(ncells-1) * multiplicity+1:degree + (ncells-2) * multiplicity + 1] - period - out[len(out) - degree + multiplicity-1:] = out[degree+multiplicity:2*degree + 1] + period + # else: + for i in range(1, ncells): + out[degree + 1 + (i - 1) * multiplicity:degree + 1 + i * multiplicity] = breaks[i] + out[0:degree + 1] = breaks[0] out[len(out) - degree - 1:] = breaks[-1] diff --git a/psydac/feec/pull_push.py b/psydac/feec/pull_push.py index c73d9394c..df9cac6be 100644 --- a/psydac/feec/pull_push.py +++ b/psydac/feec/pull_push.py @@ -340,9 +340,7 @@ def push_2d_hcurl(f1, f2, eta1, eta2, F): # # Assume that f is a list/tuple of callable functions # f1, f2 = f eta = eta1, eta2 - J_inv_value = F.jacobian_inv(*eta) - f1_value = f1(*eta) f2_value = f2(*eta) diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 1bbe49ba1..d03b49a4b 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -87,7 +87,7 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul if knots is None: if multiplicity is None:multiplicity = 1 knots = make_knots( grid, degree, periodic, multiplicity ) - + if grid is None: grid = breakpoints(knots, degree) @@ -174,12 +174,14 @@ def init_interpolation( self, dtype=float ): Greville points. """ + imat = collocation_matrix( knots = self.knots, degree = self.degree, periodic = self.periodic, normalization = self.basis, - xgrid = self.greville + xgrid = self.greville, + multiplicity = self.multiplicity ) if self.periodic: @@ -208,6 +210,7 @@ def init_histopolation( self, dtype=float): the cells defined by the extended Greville points. """ + imat = histopolation_matrix( knots = self.knots, degree = self.degree, @@ -216,6 +219,8 @@ def init_histopolation( self, dtype=float): xgrid = self.ext_greville, multiplicity = self._multiplicity ) + + self.hmat= imat if self.periodic: # Convert to CSC format and compute sparse LU decomposition diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index e2e363e2d..63ee89d50 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -191,11 +191,10 @@ 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 ): - knots = space.knots degree = space.degree - span = find_span( knots, degree, x ) - + multiplicity = space.multiplicity + span = find_span( knots, degree, x , multiplicity = multiplicity) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# @@ -622,7 +621,8 @@ def eval_field_gradient( self, field, *eta , weights=None): knots = space.knots degree = space.degree - span = find_span( knots, degree, x ) + multiplicity = space.multiplicity + span = find_span( knots, degree, x , multiplicity = multiplicity) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# From 926a9f66a6f1ffa0489709ec5bef8210a5c1d862 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 29 Sep 2023 17:18:06 +0200 Subject: [PATCH 07/37] cleaning sone added lines, etc... --- psydac/api/tests/test_api_feec_2d.py | 5 +++-- psydac/core/bsplines_kernels.py | 2 -- psydac/feec/pull_push.py | 2 ++ psydac/fem/splines.py | 4 +--- 4 files changed, 6 insertions(+), 7 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 96817ec7a..45428979b 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -464,7 +464,7 @@ def discrete_energies(self, e, b): push_2d_hcurl(E.fields[0], E.fields[1], x1i, x2j, F) Bz_values[i, j] = push_2d_l2(B, x1i, x2j, F) - + # Electric field, x component fig2 = plot_field_and_error(r'E^x', x, y, Ex_values, Ex_ex(0, x, y), *gridlines) fig2.show() @@ -605,9 +605,10 @@ def discrete_energies(self, e, b): Ex_values[i, j], Ey_values[i, j] = \ push_2d_hcurl(E.fields[0], E.fields[1], x1i, x2j, F) + Bz_values[i, j] = push_2d_l2(B, x1i, x2j, F) - # ... + # Error at final time error_Ex = abs(Ex_ex(t, x, y) - Ex_values).max() error_Ey = abs(Ey_ex(t, x, y) - Ey_values).max() diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index 06489f18f..f7d565c4e 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -818,8 +818,6 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] out[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] period = breaks[-1]-breaks[0] - #out[:degree] = breaks[ncells - degree:ncells] - period - #out[len(out) - degree:] = breaks[1:degree + 1] + period out[:degree + 1] = out[ncells * multiplicity : ncells * multiplicity + degree + 1] - period out[len(out) - degree :] = out[degree+1:2*degree + 1] + period diff --git a/psydac/feec/pull_push.py b/psydac/feec/pull_push.py index df9cac6be..c73d9394c 100644 --- a/psydac/feec/pull_push.py +++ b/psydac/feec/pull_push.py @@ -340,7 +340,9 @@ def push_2d_hcurl(f1, f2, eta1, eta2, F): # # Assume that f is a list/tuple of callable functions # f1, f2 = f eta = eta1, eta2 + J_inv_value = F.jacobian_inv(*eta) + f1_value = f1(*eta) f2_value = f2(*eta) diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index d03b49a4b..56ad6a827 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -87,7 +87,7 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul if knots is None: if multiplicity is None:multiplicity = 1 knots = make_knots( grid, degree, periodic, multiplicity ) - + if grid is None: grid = breakpoints(knots, degree) @@ -174,7 +174,6 @@ def init_interpolation( self, dtype=float ): Greville points. """ - imat = collocation_matrix( knots = self.knots, degree = self.degree, @@ -210,7 +209,6 @@ def init_histopolation( self, dtype=float): the cells defined by the extended Greville points. """ - imat = histopolation_matrix( knots = self.knots, degree = self.degree, From 7f87337cfce9a1b813d856958be8af3e3a9f3bd6 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 29 Sep 2023 17:20:01 +0200 Subject: [PATCH 08/37] cleaning --- psydac/api/tests/test_api_feec_2d.py | 2 +- psydac/core/bsplines_kernels.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 45428979b..1ae0e29f2 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -605,7 +605,7 @@ def discrete_energies(self, e, b): Ex_values[i, j], Ey_values[i, j] = \ push_2d_hcurl(E.fields[0], E.fields[1], x1i, x2j, F) - + Bz_values[i, j] = push_2d_l2(B, x1i, x2j, F) # ... diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index f7d565c4e..36098b80c 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -39,7 +39,6 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, multiplicity : int = 1 high = len(knots)-multiplicity-degree # Check if point is exactly on left/right boundary, or outside domain - if x <= knots[low ]: return low if x >= knots[high]: return high-1 From 424efceeff0c24e7e2c4355fd0f7b457703c03d7 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 29 Sep 2023 17:34:52 +0200 Subject: [PATCH 09/37] added periodic flag to find_span and added it to the places where the function is called --- psydac/api/fem.py | 5 +++-- psydac/api/grid.py | 15 ++++++++------- psydac/core/bsplines.py | 8 ++++---- psydac/core/bsplines_kernels.py | 14 +++++++++----- psydac/core/tests/test_bsplines.py | 14 +++++++------- psydac/core/tests/test_bsplines_pyccel.py | 8 ++++---- psydac/fem/splines.py | 2 +- psydac/fem/tensor.py | 6 ++++-- 8 files changed, 40 insertions(+), 32 deletions(-) diff --git a/psydac/api/fem.py b/psydac/api/fem.py index 098b80904..9d221540f 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -609,7 +609,8 @@ def construct_arguments(self, with_openmp=False): if axis is None:continue space = spaces[i].spaces[axis] points_i = points[i][axis] - local_span = find_span(space.knots, space.degree, points_i[0, 0]) + periodic = space.periodic + local_span = find_span(space.knots, space.degree, periodic, points_i[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points_i[0, 0], local_span, nderiv, space.basis) map_basis[i][axis] = map_basis[i][axis].copy() @@ -1213,7 +1214,7 @@ def construct_arguments(self, with_openmp=False): nderiv = self.max_nderiv space = space.spaces[axis] points = points[axis] - local_span = find_span(space.knots, space.degree, points[0, 0]) + local_span = find_span(space.knots, space.degree, space.periodic, points[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) map_basis[axis] = map_basis[axis].copy() diff --git a/psydac/api/grid.py b/psydac/api/grid.py index 6e8b3f118..6719d9ff2 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -206,7 +206,7 @@ def __init__( self, V, nderiv , trial=False, grid=None): for i,Vi in enumerate(V): space = Vi.spaces[axis] points = grid.points[axis] - local_span = find_span(space.knots, space.degree, points[0, 0]) + local_span = find_span(space.knots, space.degree, space.periodic, points[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) @@ -230,11 +230,12 @@ def space(self): # TODO have a parallel version of this function, as done for fem def create_collocation_basis( glob_points, space, nderiv=1 ): - T = space.knots # knots sequence - p = space.degree # spline degree - n = space.nbasis # total number of control points - grid = space.breaks # breakpoints - nc = space.ncells # number of cells in domain (nc=len(grid)-1) + T = space.knots # knots sequence + p = space.degree # spline degree + n = space.nbasis # total number of control points + grid = space.breaks # breakpoints + nc = space.ncells # number of cells in domain (nc=len(grid)-1) + perio = space.perio #------------------------------------------- # GLOBAL GRID @@ -246,7 +247,7 @@ def create_collocation_basis( glob_points, space, nderiv=1 ): # glob_basis = np.zeros( (p+1,nderiv+1,nq) ) # TODO use this for local basis fct glob_basis = np.zeros( (n+p,nderiv+1,nq) ) # n+p for ghosts for iq,xq in enumerate(glob_points): - span = find_span( T, p, xq ) + span = find_span( T, p, perio, xq ) glob_spans[iq] = span ders = basis_funs_all_ders( T, p, xq, span, nderiv ) diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 4803c8bfa..84770cd3b 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -57,7 +57,7 @@ #============================================================================== -def find_span(knots, degree, x, multiplicity = 1): +def find_span(knots, degree, x, periodic, multiplicity = 1): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -84,10 +84,10 @@ def find_span(knots, degree, x, multiplicity = 1): x = float(x) knots = np.ascontiguousarray(knots, dtype=float) multiplicity = int(multiplicity) - return find_span_p(knots, degree, x, multiplicity = multiplicity) + return find_span_p(knots, degree, x, periodic, multiplicity = multiplicity) #============================================================================== -def find_spans(knots, degree, x, out=None): +def find_spans(knots, degree, x, periodic, out=None): """ Determine the knot span index at a set of locations x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -122,7 +122,7 @@ def find_spans(knots, degree, x, out=None): else: assert out.shape == x.shape and out.dtype == np.dtype('int') - find_spans_p(knots, degree, x, out) + find_spans_p(knots, degree, x, periodic, out) return out #============================================================================== diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index 36098b80c..88f6f8f0c 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -5,7 +5,7 @@ import numpy as np # ============================================================================= -def find_span_p(knots: 'float[:]', degree: int, x: float, multiplicity : int = 1): +def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multiplicity : int = 1): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -36,7 +36,11 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, multiplicity : int = 1 """ # Knot index at left/right boundary low = degree - high = len(knots)-multiplicity-degree + if periodic : + high = len(knots)-multiplicity-degree + else : + high = len(knots)-1-degree + # Check if point is exactly on left/right boundary, or outside domain if x <= knots[low ]: return low @@ -55,7 +59,7 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, multiplicity : int = 1 # ============================================================================= -def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', multiplicity : int = 1): +def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', periodic : bool, multiplicity : int = 1): """ Determine the knot span index at a set of locations x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -81,7 +85,7 @@ def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', m n = x.shape[0] for i in range(n): - out[i] = find_span_p(knots, degree, x[i], multiplicity=multiplicity) + out[i] = find_span_p(knots, degree, x[i], periodic, multiplicity=multiplicity) # ============================================================================= @@ -421,7 +425,7 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali basis = np.zeros((nx, degree + 1)) spans = np.zeros(nx, dtype=int) - find_spans_p(knots, degree, xgrid, spans, multiplicity = multiplicity) + find_spans_p(knots, degree, xgrid, spans, periodic, multiplicity = multiplicity) basis_funs_array_p(knots, degree, xgrid, spans, basis) # Fill in non-zero matrix values diff --git a/psydac/core/tests/test_bsplines.py b/psydac/core/tests/test_bsplines.py index 4b3f0b987..e3bfe6666 100644 --- a/psydac/core/tests/test_bsplines.py +++ b/psydac/core/tests/test_bsplines.py @@ -32,9 +32,9 @@ def test_find_span( lims, nc, p, eps=1e-12 ): knots = np.r_[ [grid[0]]*p, grid, [grid[-1]]*p ] for i,xi in enumerate( grid ): - assert find_span( knots, p, x=xi-eps ) == p + max( 0, i-1 ) - 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 ) + assert find_span( knots, p, False, x=xi-eps ) == p + max( 0, i-1 ) + assert find_span( knots, p, False, x=xi ) == p + min( i, nc-1 ) + assert find_span( knots, p, False, x=xi+eps ) == p + min( i, nc-1 ) #============================================================================== @pytest.mark.parametrize( 'lims', ([0,1], [-2,3]) ) @@ -48,7 +48,7 @@ def test_basis_funs( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x ) + span = find_span( knots, p, False, x ) basis = basis_funs( knots, p, x, span ) assert len( basis ) == p+1 assert np.all( basis >= 0 ) @@ -66,7 +66,7 @@ def test_basis_funs_1st_der( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x ) + span = find_span( knots, p, False, x ) ders = basis_funs_1st_der( knots, p, x, span ) assert len( ders ) == p+1 assert abs( sum( ders ) ) < tol @@ -86,7 +86,7 @@ def test_basis_funs_all_ders( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x ) + span = find_span( knots, p, False, x ) ders = basis_funs_all_ders( knots, p, x, span, n ) # Test output array @@ -195,7 +195,7 @@ def test_cell_index(i_grid, expected): yy = np.zeros( (len(xx), nb) ) zz = np.zeros( (len(xx), nb) ) for i,x in enumerate( xx ): - span = find_span( knots, p, x ) + span = find_span( knots, p, False, x ) yy[i,span-p:span+1] = basis_funs ( knots, p, x, span ) zz[i,span-p:span+1] = basis_funs_1st_der( knots, p, x, span ) diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index f1bc9129b..1a6206879 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -493,7 +493,7 @@ def basis_integrals_true(knots, degree): @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_find_span(knots, degree, x): expected = find_span_true(knots, degree, x) - out = find_span(knots, degree, x) + out = find_span(knots, degree, False, x) assert np.allclose(expected, out, atol=ATOL, rtol=RTOL) @@ -511,7 +511,7 @@ def test_find_span(knots, degree, x): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs(knots, degree, x): - span = find_span(knots, degree, x) + span = find_span(knots, degree, False, x) expected = basis_funs_true(knots, degree, x, span) out = basis_funs(knots, degree, x, span) @@ -531,7 +531,7 @@ def test_basis_funs(knots, degree, x): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs_1st_der(knots, degree, x): - span = find_span(knots, degree, x) + span = find_span(knots, degree, False, x) expected = basis_funs_1st_der_true(knots, degree, x, span) out = basis_funs_1st_der(knots, degree, x, span) @@ -553,7 +553,7 @@ def test_basis_funs_1st_der(knots, degree, x): @pytest.mark.parametrize('n', (2, 3, 4, 5)) @pytest.mark.parametrize('normalization', ('B', 'M')) def test_basis_funs_all_ders(knots, degree, x, n, normalization): - span = find_span(knots, degree, x) + span = find_span(knots, degree, False, x) expected = basis_funs_all_ders_true(knots, degree, x, span, n, normalization) out = basis_funs_all_ders(knots, degree, x, span, n, normalization) diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 56ad6a827..62a4623f7 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -292,7 +292,7 @@ def eval_field(self, field, *eta , weights=None): eta = eta[0] - span = find_span( self.knots, self.degree, eta) + span = find_span( self.knots, self.degree, self.periodic, eta) basis_array = basis_funs( self.knots, self.degree, eta, span) index = slice(span-self.degree, span + 1) diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 63ee89d50..6a9b0b86b 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -194,7 +194,8 @@ def eval_field( self, field, *eta, weights=None): knots = space.knots degree = space.degree multiplicity = space.multiplicity - span = find_span( knots, degree, x , multiplicity = multiplicity) + periodic = space.periodic + span = find_span( knots, degree, x, periodic, multiplicity = multiplicity) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# @@ -622,7 +623,8 @@ def eval_field_gradient( self, field, *eta , weights=None): knots = space.knots degree = space.degree multiplicity = space.multiplicity - span = find_span( knots, degree, x , multiplicity = multiplicity) + periodic = space.periodic + span = find_span( knots, degree, x, periodic, multiplicity = multiplicity) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# From fcf6f92011f39676dcdc339afe6dc46796efe7ea Mon Sep 17 00:00:00 2001 From: vcarlier Date: Mon, 2 Oct 2023 19:13:13 +0200 Subject: [PATCH 10/37] solving a few problem in the tests --- psydac/core/bsplines.py | 7 +++++++ psydac/core/tests/test_bsplines.py | 12 ++++++------ psydac/core/tests/test_bsplines_pyccel.py | 8 ++++---- 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 84770cd3b..036027fc0 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -75,6 +75,12 @@ def find_span(knots, degree, x, periodic, multiplicity = 1): x : float Location of interest. + + periodic : bool + Indicating if the domain is periodic. + + multiplicity : int + Multplicity in the knot sequence. Returns ------- @@ -84,6 +90,7 @@ def find_span(knots, degree, x, periodic, multiplicity = 1): x = float(x) knots = np.ascontiguousarray(knots, dtype=float) multiplicity = int(multiplicity) + periodic = bool(periodic) return find_span_p(knots, degree, x, periodic, multiplicity = multiplicity) #============================================================================== diff --git a/psydac/core/tests/test_bsplines.py b/psydac/core/tests/test_bsplines.py index e3bfe6666..0187fde5b 100644 --- a/psydac/core/tests/test_bsplines.py +++ b/psydac/core/tests/test_bsplines.py @@ -32,9 +32,9 @@ def test_find_span( lims, nc, p, eps=1e-12 ): knots = np.r_[ [grid[0]]*p, grid, [grid[-1]]*p ] for i,xi in enumerate( grid ): - assert find_span( knots, p, False, x=xi-eps ) == p + max( 0, i-1 ) - assert find_span( knots, p, False, x=xi ) == p + min( i, nc-1 ) - assert find_span( knots, p, False, x=xi+eps ) == p + min( i, nc-1 ) + assert find_span( knots, p, x=xi-eps, periodic=False ) == p + max( 0, i-1 ) + assert find_span( knots, p, x=xi, periodic=False ) == p + min( i, nc-1 ) + assert find_span( knots, p, x=xi+eps, periodic=False ) == p + min( i, nc-1 ) #============================================================================== @pytest.mark.parametrize( 'lims', ([0,1], [-2,3]) ) @@ -48,7 +48,7 @@ def test_basis_funs( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, False, x ) + span = find_span( knots, p, x, False ) basis = basis_funs( knots, p, x, span ) assert len( basis ) == p+1 assert np.all( basis >= 0 ) @@ -66,7 +66,7 @@ def test_basis_funs_1st_der( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, False, x ) + span = find_span( knots, p, x, False ) ders = basis_funs_1st_der( knots, p, x, span ) assert len( ders ) == p+1 assert abs( sum( ders ) ) < tol @@ -86,7 +86,7 @@ def test_basis_funs_all_ders( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, False, x ) + span = find_span( knots, p, x, False ) ders = basis_funs_all_ders( knots, p, x, span, n ) # Test output array diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index 1a6206879..ffe3fa146 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -493,7 +493,7 @@ def basis_integrals_true(knots, degree): @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_find_span(knots, degree, x): expected = find_span_true(knots, degree, x) - out = find_span(knots, degree, False, x) + out = find_span(knots, degree, x, False) assert np.allclose(expected, out, atol=ATOL, rtol=RTOL) @@ -511,7 +511,7 @@ def test_find_span(knots, degree, x): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs(knots, degree, x): - span = find_span(knots, degree, False, x) + span = find_span(knots, degree, x, False) expected = basis_funs_true(knots, degree, x, span) out = basis_funs(knots, degree, x, span) @@ -531,7 +531,7 @@ def test_basis_funs(knots, degree, x): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs_1st_der(knots, degree, x): - span = find_span(knots, degree, False, x) + span = find_span(knots, degree, x, False) expected = basis_funs_1st_der_true(knots, degree, x, span) out = basis_funs_1st_der(knots, degree, x, span) @@ -553,7 +553,7 @@ def test_basis_funs_1st_der(knots, degree, x): @pytest.mark.parametrize('n', (2, 3, 4, 5)) @pytest.mark.parametrize('normalization', ('B', 'M')) def test_basis_funs_all_ders(knots, degree, x, n, normalization): - span = find_span(knots, degree, False, x) + span = find_span(knots, degree, x, False) expected = basis_funs_all_ders_true(knots, degree, x, span, n, normalization) out = basis_funs_all_ders(knots, degree, x, span, n, normalization) From 98a27ba12fa80368616afee05fe51e3efabae1a7 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Tue, 3 Oct 2023 09:44:12 +0200 Subject: [PATCH 11/37] had switched the arguments of findspan in a lot of places --- psydac/api/fem.py | 4 ++-- psydac/api/grid.py | 4 ++-- psydac/core/tests/test_bsplines.py | 2 +- psydac/fem/splines.py | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/psydac/api/fem.py b/psydac/api/fem.py index 9d221540f..1bec7abbe 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -610,7 +610,7 @@ def construct_arguments(self, with_openmp=False): space = spaces[i].spaces[axis] points_i = points[i][axis] periodic = space.periodic - local_span = find_span(space.knots, space.degree, periodic, points_i[0, 0]) + local_span = find_span(space.knots, space.degree, points_i[0, 0], periodic) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points_i[0, 0], local_span, nderiv, space.basis) map_basis[i][axis] = map_basis[i][axis].copy() @@ -1214,7 +1214,7 @@ def construct_arguments(self, with_openmp=False): nderiv = self.max_nderiv space = space.spaces[axis] points = points[axis] - local_span = find_span(space.knots, space.degree, space.periodic, points[0, 0]) + local_span = find_span(space.knots, space.degree, points[0, 0], space.periodic) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) map_basis[axis] = map_basis[axis].copy() diff --git a/psydac/api/grid.py b/psydac/api/grid.py index 6719d9ff2..e78f41c82 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -206,7 +206,7 @@ def __init__( self, V, nderiv , trial=False, grid=None): for i,Vi in enumerate(V): space = Vi.spaces[axis] points = grid.points[axis] - local_span = find_span(space.knots, space.degree, space.periodic, points[0, 0]) + local_span = find_span(space.knots, space.degree, points[0, 0], space.periodic) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) @@ -247,7 +247,7 @@ def create_collocation_basis( glob_points, space, nderiv=1 ): # glob_basis = np.zeros( (p+1,nderiv+1,nq) ) # TODO use this for local basis fct glob_basis = np.zeros( (n+p,nderiv+1,nq) ) # n+p for ghosts for iq,xq in enumerate(glob_points): - span = find_span( T, p, perio, xq ) + span = find_span( T, p, xq, perio ) glob_spans[iq] = span ders = basis_funs_all_ders( T, p, xq, span, nderiv ) diff --git a/psydac/core/tests/test_bsplines.py b/psydac/core/tests/test_bsplines.py index 0187fde5b..c539bbd9a 100644 --- a/psydac/core/tests/test_bsplines.py +++ b/psydac/core/tests/test_bsplines.py @@ -195,7 +195,7 @@ def test_cell_index(i_grid, expected): yy = np.zeros( (len(xx), nb) ) zz = np.zeros( (len(xx), nb) ) for i,x in enumerate( xx ): - span = find_span( knots, p, False, x ) + span = find_span( knots, p, x, False ) yy[i,span-p:span+1] = basis_funs ( knots, p, x, span ) zz[i,span-p:span+1] = basis_funs_1st_der( knots, p, x, span ) diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 62a4623f7..4cdea58c4 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -292,7 +292,7 @@ def eval_field(self, field, *eta , weights=None): eta = eta[0] - span = find_span( self.knots, self.degree, self.periodic, eta) + span = find_span( self.knots, self.degree, eta, self.periodic) basis_array = basis_funs( self.knots, self.degree, eta, span) index = slice(span-self.degree, span + 1) From 7cea74a595d24a0cd2a3c6f1a2b77da1b1ad24ec Mon Sep 17 00:00:00 2001 From: vcarlier Date: Tue, 3 Oct 2023 10:20:28 +0200 Subject: [PATCH 12/37] fix make_knots_true --- psydac/core/tests/test_bsplines_pyccel.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index ffe3fa146..ceb276c28 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -377,17 +377,23 @@ def make_knots_true( breaks, degree, periodic, multiplicity=1 ): assert len(breaks) > degree p = degree - T = np.zeros( multiplicity*len(breaks[1:-1])+2+2*p ) + if periodic : + T = np.zeros(multiplicity * len(breaks[1:]) + 1 + 2 * degree) + else : + T = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) - T[p+1:-p-1] = np.repeat(breaks[1:-1], multiplicity) - T[p] = breaks[ 0] - T[-p-1] = breaks[-1] + if periodic: period = breaks[-1]-breaks[0] - T[0:p] = [xi-period for xi in breaks[-p-1:-1 ]] - T[-p:] = [xi+period for xi in breaks[ 1:p+1]] + T[p+1:-p] = np.repeat(breaks[1:], multiplicity) + + T[0:p+1] = [xi-period for xi in T[-2*p-1:-p ]] + T[-p:] = [xi+period for xi in T[ p+1:2*p+1]] else: + T[p+1:-p-1] = np.repeat(breaks[1:-1], multiplicity) + T[p] = breaks[ 0] + T[-p-1] = breaks[-1] T[0:p+1] = breaks[ 0] T[-p-1:] = breaks[-1] From eae659af83eafd88e01159a3cd1f22d30a1f6e84 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Tue, 3 Oct 2023 11:46:37 +0200 Subject: [PATCH 13/37] add perio_hist flag to collocation_matrix_p to compute the right spans in the periodic case --- psydac/api/tests/test_api_feec_2d.py | 24 ++++++++++++++++++++---- psydac/core/bsplines_kernels.py | 21 ++++++++++++++++----- psydac/fem/splines.py | 1 - 3 files changed, 36 insertions(+), 10 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 1ae0e29f2..ea2a779b7 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -729,13 +729,22 @@ def test_maxwell_2d_multiplicity(): mult = 2 ) + TOL = 1e-6 + ref = dict(error_l2_Ex = 4.0298626142e-02, + error_l2_Ey = 4.0298626142e-02, + error_l2_Bz = 3.0892842064e-01) + + assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL + assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL + assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL + def test_maxwell_2d_periodic_multiplicity(): namespace = run_maxwell_2d_TE( use_spline_mapping = False, eps = 0.5, - ncells = 3, - degree = 3, + ncells = 10, + degree = 5, periodic = True, Cp = 0.5, nsteps = 1, @@ -747,6 +756,15 @@ def test_maxwell_2d_periodic_multiplicity(): verbose = False, mult =2 ) + + TOL = 1e-6 + ref = dict(error_l2_Ex = 2.0425545241e-02, + error_l2_Ey = 2.0425545241e-02, + error_l2_Bz = 1.2401168441e-01) + + assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL + assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL + assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL def test_maxwell_2d_dirichlet(): @@ -866,8 +884,6 @@ def test_maxwell_2d_dirichlet_par(): #============================================================================== if __name__ == '__main__': - test_maxwell_2d_periodic_multiplicity() - exit() import argparse parser = argparse.ArgumentParser( diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index 88f6f8f0c..1cbc2b049 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -23,7 +23,13 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi x : float Location of interest. - + + periodic : bool + Indicating if the domain is periodic. + + multiplicity : int + Multplicity in the knot sequence. + Returns ------- span : int @@ -40,7 +46,6 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi high = len(knots)-multiplicity-degree else : high = len(knots)-1-degree - # Check if point is exactly on left/right boundary, or outside domain if x <= knots[low ]: return low @@ -387,7 +392,7 @@ def basis_integrals_p(knots: 'float[:]', degree: int, out: 'float[:]'): # ============================================================================= def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normalization: bool, xgrid: 'float[:]', - out: 'float[:,:]', multiplicity : int = 1): + out: 'float[:,:]', multiplicity : int = 1, perio_hist : bool = False): """ Compute the collocation matrix :math:`C_ij = B_j(x_i)`, which contains the values of each B-spline basis function :math:`B_j` at all locations :math:`x_i`. @@ -414,6 +419,9 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali out : array The result will be inserted into this array. It should be of the appropriate shape and dtype. + + perio_hist : bool + True if this function was called by histopolation_matrix_p with a periodic domain. """ # Number of basis functions (in periodic case remove degree repeated elements) nb = len(knots)-degree-1 @@ -425,7 +433,9 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali basis = np.zeros((nx, degree + 1)) spans = np.zeros(nx, dtype=int) - find_spans_p(knots, degree, xgrid, spans, periodic, multiplicity = multiplicity) + #use perio_hist to have the right spans in the case this function was called by histopolation_matrix_p + #with a periodic domain (else find_spans will return wrong results for points at the boundary) + find_spans_p(knots, degree, xgrid, spans, (periodic or perio_hist), multiplicity = multiplicity) basis_funs_array_p(knots, degree, xgrid, spans, basis) # Fill in non-zero matrix values @@ -549,7 +559,8 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma False, xgrid_new[:actual_len], colloc, - multiplicity = multiplicity) + multiplicity = multiplicity, + perio_hist = periodic) m = colloc.shape[0] - 1 n = colloc.shape[1] - 1 diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 4cdea58c4..1b3e1a898 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -218,7 +218,6 @@ def init_histopolation( self, dtype=float): multiplicity = self._multiplicity ) - self.hmat= imat if self.periodic: # Convert to CSC format and compute sparse LU decomposition From f4242f4a508645b2fcb9983bc926e3c15e359596 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Tue, 3 Oct 2023 11:49:20 +0200 Subject: [PATCH 14/37] forget periodic flag in _refinement_matrix_one_stage --- psydac/core/bsplines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 036027fc0..9d004a91f 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -1008,7 +1008,7 @@ def alpha_function(i, k, t, n, p, knots): mat = np.zeros((n+1,n)) - left = find_span( knots, p, t ) + left = find_span( knots, p, t, False ) # ... j = 0 From 90bb3168e3318693507646e2ca85d387f5dab47c Mon Sep 17 00:00:00 2001 From: vcarlier Date: Tue, 3 Oct 2023 14:01:09 +0200 Subject: [PATCH 15/37] update the error in test --- psydac/api/tests/test_api_feec_2d.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index ea2a779b7..82ae9bb99 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -626,9 +626,9 @@ def discrete_energies(self, e, b): error_l2_Ex = np.sqrt(derham_h.V1.spaces[0].integral(errx)) error_l2_Ey = np.sqrt(derham_h.V1.spaces[1].integral(erry)) error_l2_Bz = np.sqrt(derham_h.V0.integral(errz)) - print('L2 norm of error on Ex(t,x,y) at final time: {:.2e}'.format(error_l2_Ex)) - print('L2 norm of error on Ey(t,x,y) at final time: {:.2e}'.format(error_l2_Ey)) - print('L2 norm of error on Bz(t,x,y) at final time: {:.2e}'.format(error_l2_Bz)) + print('L2 norm of error on Ex(t,x,y) at final time: {:.10e}'.format(error_l2_Ex)) + print('L2 norm of error on Ey(t,x,y) at final time: {:.10e}'.format(error_l2_Ey)) + print('L2 norm of error on Bz(t,x,y) at final time: {:.10e}'.format(error_l2_Bz)) if diagnostics_interval: @@ -758,9 +758,9 @@ def test_maxwell_2d_periodic_multiplicity(): ) TOL = 1e-6 - ref = dict(error_l2_Ex = 2.0425545241e-02, - error_l2_Ey = 2.0425545241e-02, - error_l2_Bz = 1.2401168441e-01) + ref = dict(error_l2_Ex = 8.2398553953e-02, + error_l2_Ey = 8.2398553953e-02, + error_l2_Bz = 2.6713526134e-01) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL @@ -884,6 +884,9 @@ def test_maxwell_2d_dirichlet_par(): #============================================================================== if __name__ == '__main__': + test_maxwell_2d_periodic_multiplicity() + exit() + import argparse parser = argparse.ArgumentParser( From 1f1a904eef658497678260ee96dd3f7e336ee907 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 4 Oct 2023 14:05:12 +0200 Subject: [PATCH 16/37] forgot to remove change in print --- psydac/api/tests/test_api_feec_2d.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 82ae9bb99..32b86ffd1 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -626,9 +626,9 @@ def discrete_energies(self, e, b): error_l2_Ex = np.sqrt(derham_h.V1.spaces[0].integral(errx)) error_l2_Ey = np.sqrt(derham_h.V1.spaces[1].integral(erry)) error_l2_Bz = np.sqrt(derham_h.V0.integral(errz)) - print('L2 norm of error on Ex(t,x,y) at final time: {:.10e}'.format(error_l2_Ex)) - print('L2 norm of error on Ey(t,x,y) at final time: {:.10e}'.format(error_l2_Ey)) - print('L2 norm of error on Bz(t,x,y) at final time: {:.10e}'.format(error_l2_Bz)) + print('L2 norm of error on Ex(t,x,y) at final time: {:.2e}'.format(error_l2_Ex)) + print('L2 norm of error on Ey(t,x,y) at final time: {:.2e}'.format(error_l2_Ey)) + print('L2 norm of error on Bz(t,x,y) at final time: {:.2e}'.format(error_l2_Bz)) if diagnostics_interval: From 540a7abfe51e2c07721d3b7c18168fa8be3232b2 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Tue, 10 Oct 2023 18:27:15 +0200 Subject: [PATCH 17/37] add tests in 1d and 3d, found a problem probably comming from creation of greville pts --- psydac/api/tests/test_api_feec_1d.py | 31 +++++++++- psydac/api/tests/test_api_feec_3d.py | 93 ++++++++++++++++++++++++++-- 2 files changed, 118 insertions(+), 6 deletions(-) diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index b94013cac..dbac9b768 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -46,7 +46,7 @@ def update_plot(ax, t, sol_ex, sol_num): #============================================================================== def run_maxwell_1d(*, L, eps, ncells, degree, periodic, Cp, nsteps, tend, splitting_order, plot_interval, diagnostics_interval, - bc_mode, tol, verbose): + bc_mode, tol, verbose, mult=1): import numpy as np import matplotlib.pyplot as plt @@ -132,7 +132,7 @@ class CollelaMapping1D(Mapping): # Discrete physical domain and discrete DeRham sequence domain_h = discretize(domain, ncells=[ncells], periodic=[periodic], comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=[degree]) + derham_h = discretize(derham, domain_h, degree=[degree], multiplicity=[mult]) # Discrete bilinear forms a0_h = discretize(a0, domain_h, (derham_h.V0, derham_h.V0), backend=PSYDAC_BACKEND_PYTHON) @@ -496,6 +496,33 @@ def test_maxwell_1d_periodic(): assert abs(namespace['error_E'] - ref['error_E']) / ref['error_E'] <= TOL assert abs(namespace['error_B'] - ref['error_B']) / ref['error_B'] <= TOL + +def test_maxwell_1d_periodic_mult(): + + namespace = run_maxwell_1d( + L = 1.0, + eps = 0.5, + ncells = 30, + degree = 3, + periodic = True, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + bc_mode = None, + verbose = False, + mult = 2 + ) + + TOL = 1e-6 + ref = dict(error_E = 4.24689277e-04, + error_B = 4.03195666e-04) + + assert abs(namespace['error_E'] - ref['error_E']) / ref['error_E'] <= TOL + assert abs(namespace['error_B'] - ref['error_B']) / ref['error_B'] <= TOL def test_maxwell_1d_dirichlet_strong(): diff --git a/psydac/api/tests/test_api_feec_3d.py b/psydac/api/tests/test_api_feec_3d.py index 323d23c41..0da8a8b76 100644 --- a/psydac/api/tests/test_api_feec_3d.py +++ b/psydac/api/tests/test_api_feec_3d.py @@ -86,7 +86,7 @@ def evaluation_all_times(fields, x, y, z): return ak_value #================================================================================== -def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter): +def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=1): #------------------------------------------------------------------------------ # Symbolic objects: SymPDE @@ -113,7 +113,7 @@ def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, pe #------------------------------------------------------------------------------ domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=degree) + derham_h = discretize(derham, domain_h, degree=degree, multiplicity = [mult,mult,mult]) a1_h = discretize(a1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL) a2_h = discretize(a2, domain_h, (derham_h.V2, derham_h.V2), backend=PSYDAC_BACKEND_GPYCCEL) @@ -181,7 +181,7 @@ def run_maxwell_3d_scipy(logical_domain, mapping, e_ex, b_ex, ncells, degree, pe return error #================================================================================== -def run_maxwell_3d_stencil(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter): +def run_maxwell_3d_stencil(logical_domain, mapping, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=1): #------------------------------------------------------------------------------ # Symbolic objects: SymPDE @@ -208,7 +208,7 @@ def run_maxwell_3d_stencil(logical_domain, mapping, e_ex, b_ex, ncells, degree, #------------------------------------------------------------------------------ domain_h = discretize(domain, ncells=ncells, periodic=periodic, comm=MPI.COMM_WORLD) - derham_h = discretize(derham, domain_h, degree=degree) + derham_h = discretize(derham, domain_h, degree=degree, multiplicity = [mult,mult,mult]) a1_h = discretize(a1, domain_h, (derham_h.V1, derham_h.V1), backend=PSYDAC_BACKEND_GPYCCEL) a2_h = discretize(a2, domain_h, (derham_h.V2, derham_h.V2), backend=PSYDAC_BACKEND_GPYCCEL) @@ -317,6 +317,46 @@ class CollelaMapping3D(Mapping): error = run_maxwell_3d_scipy(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter) assert abs(error - 0.04294761712765949) < 1e-9 +#------------------------------------------------------------------------------ +def test_maxwell_3d_1_mult(): + class CollelaMapping3D(Mapping): + + _expressions = {'x': 'k1*(x1 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', + 'y': 'k2*(x2 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', + 'z': 'k3*x3'} + + _ldim = 3 + _pdim = 3 + + M = CollelaMapping3D('M', k1=1, k2=1, k3=1, eps=0.1) + logical_domain = Cube('C', bounds1=(0, 1), bounds2=(0, 1), bounds3=(0, 1)) + + # exact solution + e_ex_0 = lambda t, x, y, z: 0 + e_ex_1 = lambda t, x, y, z: -np.cos(2*np.pi*t-2*np.pi*z) + e_ex_2 = lambda t, x, y, z: 0 + + e_ex = (e_ex_0, e_ex_1, e_ex_2) + + b_ex_0 = lambda t, x, y, z : np.cos(2*np.pi*t-2*np.pi*z) + b_ex_1 = lambda t, x, y, z : 0 + b_ex_2 = lambda t, x, y, z : 0 + + b_ex = (b_ex_0, b_ex_1, b_ex_2) + + #space parameters + ncells = [2**4, 2**3, 2**5] + degree = [2, 2, 2] + periodic = [True, True, True] + + #time parameters + dt = 0.5*1/max(ncells) + niter = 10 + T = dt*niter + + error = run_maxwell_3d_scipy(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult = 2) + assert abs(error - 0.04294761712765949) < 1e-9 + #------------------------------------------------------------------------------ def test_maxwell_3d_2(): class CollelaMapping3D(Mapping): @@ -356,6 +396,46 @@ class CollelaMapping3D(Mapping): error = run_maxwell_3d_stencil(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter) assert abs(error - 0.24586986658559362) < 1e-9 + +#------------------------------------------------------------------------------ +def test_maxwell_3d_2_mult(): + class CollelaMapping3D(Mapping): + + _expressions = {'x': 'k1*(x1 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', + 'y': 'k2*(x2 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', + 'z': 'k3*x3'} + + _ldim = 3 + _pdim = 3 + + M = CollelaMapping3D('M', k1=1, k2=1, k3=1, eps=0.1) + logical_domain = Cube('C', bounds1=(0, 1), bounds2=(0, 1), bounds3=(0, 1)) + + # exact solution + e_ex_0 = lambda t, x, y, z: 0 + e_ex_1 = lambda t, x, y, z: -np.cos(2*np.pi*t-2*np.pi*z) + e_ex_2 = lambda t, x, y, z: 0 + + e_ex = (e_ex_0, e_ex_1, e_ex_2) + + b_ex_0 = lambda t, x, y, z : np.cos(2*np.pi*t-2*np.pi*z) + b_ex_1 = lambda t, x, y, z : 0 + b_ex_2 = lambda t, x, y, z : 0 + + b_ex = (b_ex_0, b_ex_1, b_ex_2) + + #space parameters + ncells = [7, 7, 7] + degree = [2, 2, 2] + periodic = [True, True, True] + + #time parameters + dt = 0.5*1/max(ncells) + niter = 2 + T = dt*niter + + error = run_maxwell_3d_stencil(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=2) + assert abs(error - 0.24586986658559362) < 1e-9 #============================================================================== # CLEAN UP SYMPY NAMESPACE @@ -368,4 +448,9 @@ def teardown_module(): def teardown_function(): from sympy.core import cache cache.clear_cache() + +if __name__ == '__main__' : + test_maxwell_3d_1_mult() + test_maxwell_3d_2_mult() + From bcc04db2d7ea8e92cb66c3c55107cbed831d3b83 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 11 Oct 2023 13:19:39 +0200 Subject: [PATCH 18/37] changes in greville in case of minimal regularity --- psydac/api/tests/test_api_feec_2d.py | 31 +++++++++++++++++++++++++++- psydac/core/bsplines.py | 5 +++-- psydac/core/bsplines_kernels.py | 6 +++--- psydac/fem/splines.py | 4 ++-- 4 files changed, 38 insertions(+), 8 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 32b86ffd1..6d2f983fd 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -765,6 +765,35 @@ def test_maxwell_2d_periodic_multiplicity(): assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL + + +def test_maxwell_2d_periodic_multiplicity_equal_deg(): + + namespace = run_maxwell_2d_TE( + use_spline_mapping = False, + eps = 0.5, + ncells = 10, + degree = 2, + periodic = True, + Cp = 0.5, + nsteps = 1, + tend = None, + splitting_order = 2, + plot_interval = 0, + diagnostics_interval = 0, + tol = 1e-6, + verbose = False, + mult =2 + ) + + TOL = 1e-6 + ref = dict(error_l2_Ex = 7.43229667e-02, + error_l2_Ey = 7.43229667e-02, + error_l2_Bz = 2.59863979e-01) + + assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL + assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL + assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL def test_maxwell_2d_dirichlet(): @@ -884,7 +913,7 @@ def test_maxwell_2d_dirichlet_par(): #============================================================================== if __name__ == '__main__': - test_maxwell_2d_periodic_multiplicity() + test_maxwell_2d_periodic_multiplicity_equal_deg() exit() import argparse diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 9d004a91f..1ec49209d 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -481,7 +481,7 @@ def breakpoints(knots, degree, tol=1e-15, out=None): return out[:i_final] #============================================================================== -def greville(knots, degree, periodic, out=None): +def greville(knots, degree, periodic, out=None, multiplicity=1): """ Compute coordinates of all Greville points. @@ -510,7 +510,8 @@ def greville(knots, degree, periodic, out=None): if out is None: n = len(knots) - 2 * degree - 1 if periodic else len(knots) - degree - 1 out = np.zeros(n) - greville_p(knots, degree, periodic, out) + multiplicity = int(multiplicity) + greville_p(knots, degree, periodic, out, multiplicity) return out #=============================================================================== diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index 1cbc2b049..0e7e4c45b 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -702,7 +702,7 @@ def breakpoints_p(knots: 'float[:]', degree: int, out: 'float[:]', tol: float = # ============================================================================= -def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]'): +def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]', multiplicity: int=1): """ Compute coordinates of all Greville points. @@ -726,9 +726,9 @@ def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]'): n = len(T)-2*p-1 if periodic else len(T)-p-1 # Compute greville abscissas as average of p consecutive knot values - if p == 0: + if p == multiplicity-1: for i in range(n): - out[i] = sum(T[i:i + 2]) / 2 + out[i] = sum(T[i:i + p + 2]) / (p + 2) else: for i in range(1, 1+n): out[i - 1] = sum(T[i:i + p]) / p diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 1b3e1a898..190caa50c 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -131,8 +131,8 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul self._nbasis = nbasis self._breaks = grid self._ncells = len(grid) - 1 - self._greville = greville(knots, degree, periodic) - self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic) + self._greville = greville(knots, degree, periodic, multiplicity = multiplicity) + self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic, multiplicity = multiplicity) self._scaling_array = scaling_array From 82e65048ba1a65e4b0fe78a2631fc32fe3fc92a1 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 11 Oct 2023 13:39:36 +0200 Subject: [PATCH 19/37] remove added line in test --- psydac/api/tests/test_api_feec_2d.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 6d2f983fd..822cdc3be 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -912,9 +912,6 @@ def test_maxwell_2d_dirichlet_par(): # SCRIPT CAPABILITIES #============================================================================== if __name__ == '__main__': - - test_maxwell_2d_periodic_multiplicity_equal_deg() - exit() import argparse From b03626cbc0dc25fb600e1e63acd64b970a235619 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 11 Oct 2023 14:06:36 +0200 Subject: [PATCH 20/37] add test on new find_span --- psydac/api/tests/test_api_feec_2d.py | 2 +- psydac/core/tests/test_bsplines_pyccel.py | 20 ++++++++++++-------- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 822cdc3be..4c606c14a 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -912,7 +912,7 @@ def test_maxwell_2d_dirichlet_par(): # SCRIPT CAPABILITIES #============================================================================== if __name__ == '__main__': - + import argparse parser = argparse.ArgumentParser( diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index ceb276c28..aeb9f68cd 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -21,11 +21,13 @@ # "True" Functions ############################################################################### -def find_span_true( knots, degree, x ): +def find_span_true( knots, degree, x, periodic, multiplicity=1 ): # Knot index at left/right boundary low = degree - high = len(knots)-1-degree - + if periodic : + high = len(knots)-multiplicity-degree + else : + high = len(knots)-1-degree # Check if point is exactly on left/right boundary, or outside domain if x <= knots[low ]: return low if x >= knots[high]: return high-1 @@ -209,7 +211,7 @@ def collocation_matrix_true(knots, degree, periodic, normalization, xgrid): # Fill in non-zero matrix values for i,x in enumerate( xgrid ): - span = find_span_true( knots, degree, x ) + span = find_span_true( knots, degree, x, periodic ) basis = basis_funs_true( knots, degree, x, span ) mat[i,js(span)] = normalize(basis, span) @@ -459,7 +461,7 @@ def basis_ders_on_quad_grid_true(knots, degree, quad_grid, nders, normalization) for ie in range(ne): xx = quad_grid[ie, :] for iq, xq in enumerate(xx): - span = find_span_true(knots, degree, xq) + span = find_span_true(knots, degree, xq, False) ders = basis_funs_all_ders_true(knots, degree, xq, span, nders) if normalization == 'M': ders *= scaling[None, span-degree:span+1] @@ -497,9 +499,11 @@ def basis_integrals_true(knots, degree): (np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0]), 2), (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) -def test_find_span(knots, degree, x): - expected = find_span_true(knots, degree, x) - out = find_span(knots, degree, x, False) +@pytest.mark.parametrize('periodic', [True, False]) +@pytest.mark.parametrize('multiplicity', [1, 2, 3]) +def test_find_span(knots, degree, x, periodic, multiplicity): + expected = find_span_true(knots, degree, x, periodic, multiplicity=multiplicity) + out = find_span(knots, degree, x, periodic, multiplicity=multiplicity) assert np.allclose(expected, out, atol=ATOL, rtol=RTOL) From 8dbc52fbf5b939e4ad3d2b539d9889bdecd0a114 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 11 Oct 2023 17:48:00 +0200 Subject: [PATCH 21/37] fix 3d test --- psydac/api/tests/test_api_feec_3d.py | 43 +--------------------------- 1 file changed, 1 insertion(+), 42 deletions(-) diff --git a/psydac/api/tests/test_api_feec_3d.py b/psydac/api/tests/test_api_feec_3d.py index 0da8a8b76..1bbf53318 100644 --- a/psydac/api/tests/test_api_feec_3d.py +++ b/psydac/api/tests/test_api_feec_3d.py @@ -317,46 +317,6 @@ class CollelaMapping3D(Mapping): error = run_maxwell_3d_scipy(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter) assert abs(error - 0.04294761712765949) < 1e-9 -#------------------------------------------------------------------------------ -def test_maxwell_3d_1_mult(): - class CollelaMapping3D(Mapping): - - _expressions = {'x': 'k1*(x1 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', - 'y': 'k2*(x2 + eps*sin(2.*pi*x1)*sin(2.*pi*x2))', - 'z': 'k3*x3'} - - _ldim = 3 - _pdim = 3 - - M = CollelaMapping3D('M', k1=1, k2=1, k3=1, eps=0.1) - logical_domain = Cube('C', bounds1=(0, 1), bounds2=(0, 1), bounds3=(0, 1)) - - # exact solution - e_ex_0 = lambda t, x, y, z: 0 - e_ex_1 = lambda t, x, y, z: -np.cos(2*np.pi*t-2*np.pi*z) - e_ex_2 = lambda t, x, y, z: 0 - - e_ex = (e_ex_0, e_ex_1, e_ex_2) - - b_ex_0 = lambda t, x, y, z : np.cos(2*np.pi*t-2*np.pi*z) - b_ex_1 = lambda t, x, y, z : 0 - b_ex_2 = lambda t, x, y, z : 0 - - b_ex = (b_ex_0, b_ex_1, b_ex_2) - - #space parameters - ncells = [2**4, 2**3, 2**5] - degree = [2, 2, 2] - periodic = [True, True, True] - - #time parameters - dt = 0.5*1/max(ncells) - niter = 10 - T = dt*niter - - error = run_maxwell_3d_scipy(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult = 2) - assert abs(error - 0.04294761712765949) < 1e-9 - #------------------------------------------------------------------------------ def test_maxwell_3d_2(): class CollelaMapping3D(Mapping): @@ -435,7 +395,7 @@ class CollelaMapping3D(Mapping): T = dt*niter error = run_maxwell_3d_stencil(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=2) - assert abs(error - 0.24586986658559362) < 1e-9 + assert abs(error - 0.9020674104318802) < 1e-9 #============================================================================== # CLEAN UP SYMPY NAMESPACE @@ -450,7 +410,6 @@ def teardown_function(): cache.clear_cache() if __name__ == '__main__' : - test_maxwell_3d_1_mult() test_maxwell_3d_2_mult() From de64147de40a5aaf41b2020760cd864c3821fad5 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Thu, 12 Oct 2023 10:48:51 +0200 Subject: [PATCH 22/37] docstring in find_span --- psydac/core/bsplines.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 1ec49209d..aebf9f57d 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -80,7 +80,9 @@ def find_span(knots, degree, x, periodic, multiplicity = 1): Indicating if the domain is periodic. multiplicity : int - Multplicity in the knot sequence. + Multplicity in the knot sequence. One assume that all the knots have + multiplicity multiplicity, and in the none periodic case that the + boundary knots have multiplicity p+1 (has created by make_knots) Returns ------- From eef6258ae4fb799d0b427f2de50e6a8a35a0ac8b Mon Sep 17 00:00:00 2001 From: vcarlier Date: Thu, 12 Oct 2023 14:07:59 +0200 Subject: [PATCH 23/37] add docs --- psydac/api/grid.py | 5 +++-- psydac/core/bsplines.py | 10 +++++++++- psydac/core/bsplines_kernels.py | 17 +++++++++++++++-- 3 files changed, 27 insertions(+), 5 deletions(-) diff --git a/psydac/api/grid.py b/psydac/api/grid.py index e78f41c82..1771bc131 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -229,13 +229,14 @@ def space(self): #============================================================================== # TODO have a parallel version of this function, as done for fem def create_collocation_basis( glob_points, space, nderiv=1 ): - + + #V.C 12/10/23 : Is this function really used somewhere? What is it suppose to do? T = space.knots # knots sequence p = space.degree # spline degree n = space.nbasis # total number of control points grid = space.breaks # breakpoints nc = space.ncells # number of cells in domain (nc=len(grid)-1) - perio = space.perio + perio = space.periodic # periodicity of the space #------------------------------------------- # GLOBAL GRID diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index aebf9f57d..9d76ae526 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -574,7 +574,15 @@ def make_knots(breaks, degree, periodic, multiplicity=1, out=None): Create spline knots from breakpoints, with appropriate boundary conditions. Let p be spline degree. If domain is periodic, knot sequence is extended by periodicity so that first p basis functions are identical to last p. - Otherwise, knot sequence is clamped (i.e. endpoints are repeated p times). + + This is done by first creating the "true" knot sequence (the one that we + would have if the space was really periodic), by ignoring the first + breakpoint (which is equal to the last one with periodicity). Then we add + 2p+1 more knots so that the total number of basis function is n_basis+p. + p+1 of this knots are on the left by shifting the last p+1 knots previously + created and p of them are added on the left. + + Otherwise, knot sequence is clamped (i.e. endpoints have multiplicity p+1). Parameters ---------- diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index 0e7e4c45b..f6bedcc81 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -28,7 +28,9 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi Indicating if the domain is periodic. multiplicity : int - Multplicity in the knot sequence. + Multplicity in the knot sequence. One assume that all the knots have + multiplicity multiplicity, and in the none periodic case that the + boundary knots have multiplicity p+1 (has created by make_knots) Returns ------- @@ -41,6 +43,9 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi Springer-Verlag Berlin Heidelberg GmbH, 1997. """ # Knot index at left/right boundary + #In the periodic case the p+1 knots is the last knot equal to the left + #boundary, while the first knots being the right boundary is multiplicity+p + #before the last knot of the sequence (see make knots) low = degree if periodic : high = len(knots)-multiplicity-degree @@ -803,7 +808,15 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] Create spline knots from breakpoints, with appropriate boundary conditions. Let p be spline degree. If domain is periodic, knot sequence is extended by periodicity so that first p basis functions are identical to last p. - Otherwise, knot sequence is clamped (i.e. endpoints are repeated p times). + + This is done by first creating the "true" knot sequence (the one that we + would have if the space was really periodic), by ignoring the first + breakpoint (which is equal to the last one with periodicity). Then we add + 2p+1 more knots so that the total number of basis function is n_basis+p. + p+1 of this knots are on the left by shifting the last p+1 knots previously + created and p of them are added on the left. + + Otherwise, knot sequence is clamped (i.e. endpoints have multiplicity p+1). Parameters ---------- From 94c2d8dd723c02533b223d82d251b7a324cd6dfc Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 13 Oct 2023 11:02:06 +0200 Subject: [PATCH 24/37] changing make_knots and adapting the rest --- psydac/api/tests/test_api_feec_1d.py | 8 ++--- psydac/api/tests/test_api_feec_2d.py | 10 ++++-- psydac/core/bsplines.py | 16 +++------ psydac/core/bsplines_kernels.py | 50 +++++++++++----------------- 4 files changed, 34 insertions(+), 50 deletions(-) diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index dbac9b768..93f69a45d 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -406,8 +406,8 @@ def discrete_energies(self, e, b): error_E = max(abs(E_ex(t, x) - E_values)) error_B = max(abs(B_ex(t, x) - B_values)) print() - print('Max-norm of error on E(t,x) at final time: {:.2e}'.format(error_E)) - print('Max-norm of error on B(t,x) at final time: {:.2e}'.format(error_B)) + print('Max-norm of error on E(t,x) at final time: {:.8e}'.format(error_E)) + print('Max-norm of error on B(t,x) at final time: {:.8e}'.format(error_B)) # compute L2 error as well F = mapping.get_callable_mapping() @@ -518,8 +518,8 @@ def test_maxwell_1d_periodic_mult(): ) TOL = 1e-6 - ref = dict(error_E = 4.24689277e-04, - error_B = 4.03195666e-04) + ref = dict(error_E = 4.54259221e-04, + error_B = 3.71519383e-04) assert abs(namespace['error_E'] - ref['error_E']) / ref['error_E'] <= TOL assert abs(namespace['error_B'] - ref['error_B']) / ref['error_B'] <= TOL diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 4c606c14a..265619ddb 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -743,10 +743,11 @@ def test_maxwell_2d_periodic_multiplicity(): namespace = run_maxwell_2d_TE( use_spline_mapping = False, eps = 0.5, - ncells = 10, - degree = 5, + ncells = 30, + degree = 3, periodic = True, - Cp = 0.5, + Cp = 0.5, #it's seems that with higher multiplicity we need + #to be more carefull with the CFL (with C=5. error is ) nsteps = 1, tend = None, splitting_order = 2, @@ -912,6 +913,9 @@ def test_maxwell_2d_dirichlet_par(): # SCRIPT CAPABILITIES #============================================================================== if __name__ == '__main__': + + test_maxwell_2d_periodic_multiplicity() + exit() import argparse diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 9d76ae526..2fa85aa05 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -572,15 +572,10 @@ def elements_spans(knots, degree, out=None): def make_knots(breaks, degree, periodic, multiplicity=1, out=None): """ Create spline knots from breakpoints, with appropriate boundary conditions. - Let p be spline degree. If domain is periodic, knot sequence is extended - by periodicity so that first p basis functions are identical to last p. - This is done by first creating the "true" knot sequence (the one that we - would have if the space was really periodic), by ignoring the first - breakpoint (which is equal to the last one with periodicity). Then we add - 2p+1 more knots so that the total number of basis function is n_basis+p. - p+1 of this knots are on the left by shifting the last p+1 knots previously - created and p of them are added on the left. + If domain is periodic, knot sequence is extended by periodicity to have a + total of (n_cells-1)*mult+2p+2 knots (all break points are repeated mult + time and we add p+1-mult knots by periodicity at each side). Otherwise, knot sequence is clamped (i.e. endpoints have multiplicity p+1). @@ -627,10 +622,7 @@ def make_knots(breaks, degree, periodic, multiplicity=1, out=None): breaks = np.ascontiguousarray(breaks, dtype=float) if out is None: - if periodic : - out = np.zeros(multiplicity * len(breaks[1:]) + 1 + 2 * degree) - else : - out = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) + out = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) else: assert out.shape == (multiplicity * len(breaks[1:-1]) + 2 + 2 * degree,) \ and out.dtype == np.dtype('float') diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index f6bedcc81..d5c60be6a 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -47,10 +47,7 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi #boundary, while the first knots being the right boundary is multiplicity+p #before the last knot of the sequence (see make knots) low = degree - if periodic : - high = len(knots)-multiplicity-degree - else : - high = len(knots)-1-degree + high = len(knots)-1-degree # Check if point is exactly on left/right boundary, or outside domain if x <= knots[low ]: return low @@ -397,7 +394,7 @@ def basis_integrals_p(knots: 'float[:]', degree: int, out: 'float[:]'): # ============================================================================= def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normalization: bool, xgrid: 'float[:]', - out: 'float[:,:]', multiplicity : int = 1, perio_hist : bool = False): + out: 'float[:,:]', multiplicity : int = 1): """ Compute the collocation matrix :math:`C_ij = B_j(x_i)`, which contains the values of each B-spline basis function :math:`B_j` at all locations :math:`x_i`. @@ -425,8 +422,6 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali The result will be inserted into this array. It should be of the appropriate shape and dtype. - perio_hist : bool - True if this function was called by histopolation_matrix_p with a periodic domain. """ # Number of basis functions (in periodic case remove degree repeated elements) nb = len(knots)-degree-1 @@ -438,9 +433,8 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali basis = np.zeros((nx, degree + 1)) spans = np.zeros(nx, dtype=int) - #use perio_hist to have the right spans in the case this function was called by histopolation_matrix_p - #with a periodic domain (else find_spans will return wrong results for points at the boundary) - find_spans_p(knots, degree, xgrid, spans, (periodic or perio_hist), multiplicity = multiplicity) + + find_spans_p(knots, degree, xgrid, spans, periodic, multiplicity = multiplicity) basis_funs_array_p(knots, degree, xgrid, spans, basis) # Fill in non-zero matrix values @@ -564,8 +558,7 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma False, xgrid_new[:actual_len], colloc, - multiplicity = multiplicity, - perio_hist = periodic) + multiplicity = multiplicity) m = colloc.shape[0] - 1 n = colloc.shape[1] - 1 @@ -806,15 +799,10 @@ def elements_spans_p(knots: 'float[:]', degree: int, out: 'int[:]'): def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:]', multiplicity: int = 1): """ Create spline knots from breakpoints, with appropriate boundary conditions. - Let p be spline degree. If domain is periodic, knot sequence is extended - by periodicity so that first p basis functions are identical to last p. - This is done by first creating the "true" knot sequence (the one that we - would have if the space was really periodic), by ignoring the first - breakpoint (which is equal to the last one with periodicity). Then we add - 2p+1 more knots so that the total number of basis function is n_basis+p. - p+1 of this knots are on the left by shifting the last p+1 knots previously - created and p of them are added on the left. + If domain is periodic, knot sequence is extended by periodicity to have a + total of (n_cells-1)*mult+2p+2 knots (all break points are repeated mult + time and we add p+1-mult knots by periodicity at each side). Otherwise, knot sequence is clamped (i.e. endpoints have multiplicity p+1). @@ -840,22 +828,22 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] """ ncells = len(breaks) - 1 + for i in range(0, ncells+1): + out[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] + + len_out = len(out) + if periodic: - for i in range(1, ncells+1): - out[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] period = breaks[-1]-breaks[0] - out[:degree + 1] = out[ncells * multiplicity : ncells * multiplicity + degree + 1] - period - out[len(out) - degree :] = out[degree+1:2*degree + 1] + period + out[: degree + 1 - multiplicity] = out[len_out - 2 * (degree + 1 )+ multiplicity: len_out - (degree + 1)] - period + out[len(out) - (degree + 1 - multiplicity) :] = out[degree + 1:2*(degree + 1)- multiplicity] + period # else: - for i in range(1, ncells): - out[degree + 1 + (i - 1) * multiplicity:degree + 1 + i * multiplicity] = breaks[i] - - out[0:degree + 1] = breaks[0] - out[len(out) - degree - 1:] = breaks[-1] + out[0:degree + 1 - multiplicity] = breaks[0] + out[len_out - degree - 1 + multiplicity:] = breaks[-1] # ============================================================================= @@ -897,8 +885,8 @@ def elevate_knots_p(knots: 'float[:]', degree: int, periodic: bool, out: 'float[ if periodic: T, p = knots, degree period = T[len(knots) -1 - p] - T[p] - left = T[len(knots) -2 - 2 * p] - period - right = T[2 * p + 1] + period + left = T[len(knots) -2 - 2 * p + multiplicity-1] - period + right = T[2 * p + 2 - multiplicity] + period out[0] = left out[-1] = right From 991e3e563bea8de6e5bb6eae01a4f7b5ed32eb83 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 13 Oct 2023 11:02:31 +0200 Subject: [PATCH 25/37] changing make_knots and adapting the rest --- psydac/fem/splines.py | 1 - 1 file changed, 1 deletion(-) diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 190caa50c..9a7d3d933 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -133,7 +133,6 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul self._ncells = len(grid) - 1 self._greville = greville(knots, degree, periodic, multiplicity = multiplicity) self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic, multiplicity = multiplicity) - self._scaling_array = scaling_array self._parent_multiplicity = parent_multiplicity From 78b415a51987bb74a817e645987036be6e3eedca Mon Sep 17 00:00:00 2001 From: vcarlier Date: Fri, 20 Oct 2023 08:42:36 +0200 Subject: [PATCH 26/37] new tests on feec with higher multiplicity --- psydac/api/tests/test_api_feec_2d.py | 2 +- psydac/api/tests/test_assembly.py | 58 ++++++++++++++++++- psydac/feec/derivatives.py | 9 +-- .../tests/test_differentiation_matrices.py | 41 ++++++++----- psydac/feec/tests/test_global_projectors.py | 20 ++++--- 5 files changed, 102 insertions(+), 28 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 265619ddb..c0c8f0c0d 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -914,7 +914,7 @@ def test_maxwell_2d_dirichlet_par(): #============================================================================== if __name__ == '__main__': - test_maxwell_2d_periodic_multiplicity() + test_maxwell_2d_multiplicity() exit() import argparse diff --git a/psydac/api/tests/test_assembly.py b/psydac/api/tests/test_assembly.py index 0dfae99f5..b598eb55d 100644 --- a/psydac/api/tests/test_assembly.py +++ b/psydac/api/tests/test_assembly.py @@ -85,6 +85,60 @@ def test_field_and_constant(backend, dtype): assert abs(b.toarray().sum() - res) < 1e-12 print("PASSED") +#============================================================================== +def test_field_and_constant_mult(backend, dtype): + + # If 'backend' is specified, accelerate Python code by passing **kwargs + # to discretization of bilinear forms, linear forms and functionals. + kwargs = {'backend': PSYDAC_BACKENDS[backend]} if backend else {} + + domain = Square() + V = ScalarFunctionSpace('V', domain) + + # TODO: remove codomain_type when It is implemented in sympde + u = element_of(V, name='u') + v = element_of(V, name='v') + f = element_of(V, name='f') + + if dtype == 'complex': + c = Constant(name='c', complex=True) + V.codomain_type = dtype + g = I * c * f**2 + res = 1.j + cst=complex(1.0) + else: + c = Constant(name='c', real=True) + g = c * f**2 + res = 1 + cst=1.0 + + a = BilinearForm((u, v), integral(domain, g * u * v)) + l = LinearForm(v, integral(domain, g * v)) + + ncells = (5, 5) + degree = (3, 3) + multiplicity = (2,2) + domain_h = discretize(domain, ncells=ncells) + Vh = discretize(V, domain_h, degree=degree, multiplicity = multiplicity) + ah = discretize(a, domain_h, [Vh, Vh], **kwargs) + lh = discretize(l, domain_h, Vh , **kwargs) + + fh = FemField(Vh) + fh.coeffs[:] = 1 + + # Assembly call should not crash if correct arguments are used + A = ah.assemble(c=cst, f=fh) + b = lh.assemble(f=fh, c=cst) + + # Test matrix A + x = fh.coeffs + #TODO change res into np.conj(res) when the conjugate is applied in the dot product in sympde + assert abs(x.dot(A.dot(x)) - res) < 1e-12 + + # Test vector b + assert abs(b.toarray().sum() - res) < 1e-12 + print("PASSED") + #============================================================================== def test_bilinearForm_complex(backend): @@ -532,7 +586,9 @@ def test_assembly_no_synchr_args(backend): #============================================================================== if __name__ == '__main__': - test_field_and_constant(None) + test_field_and_constant_mult('pyccel-gcc','real') + exit() + test_field_and_constant(None,'real') test_multiple_fields(None) test_math_imports(None) test_non_symmetric_BilinearForm(None) diff --git a/psydac/feec/derivatives.py b/psydac/feec/derivatives.py index 28fec0c11..94dd2643d 100644 --- a/psydac/feec/derivatives.py +++ b/psydac/feec/derivatives.py @@ -101,8 +101,8 @@ def __init__(self, V, W, diffdir, *, negative=False, transposed=False): self._codomain = W # the local area in the codomain without padding - self._idslice = tuple([slice(pad, e-s+1+pad) for pad, s, e - in zip(self._codomain.pads, self._codomain.starts, self._codomain.ends)]) + self._idslice = tuple([slice(pad*m, e-s+1+pad*m) for pad, s, e, m + in zip(self._codomain.pads, self._codomain.starts, self._codomain.ends, self._codomain.shifts)]) # prepare the slices (they are of the right size then, we checked this already) # identity slice @@ -112,12 +112,13 @@ def __init__(self, V, W, diffdir, *, negative=False, transposed=False): diff_pad = self._codomain.pads[self._diffdir] diff_s = self._codomain.starts[self._diffdir] diff_e = self._codomain.ends[self._diffdir] + diff_m = self._codomain.shifts[self._diffdir] # the diffslice depends on the transposition if self._transposed: - diff_partslice = slice(diff_pad-1, diff_e-diff_s+1+diff_pad-1) + diff_partslice = slice(diff_m*diff_pad-1, diff_e-diff_s+1+diff_pad*diff_m-1) else: - diff_partslice = slice(diff_pad+1, diff_e-diff_s+1+diff_pad+1) + diff_partslice = slice(diff_m*diff_pad+1, diff_e-diff_s+1+diff_m*diff_pad+1) diffslice = tuple([diff_partslice if i==self._diffdir else idslice[i] for i in range(self._domain.ndim)]) diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 44f7b91cc..3feccd55d 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -12,6 +12,8 @@ from psydac.feec.derivatives import ScalarCurl_2D, VectorCurl_2D, Curl_3D from psydac.feec.derivatives import Divergence_2D, Divergence_3D +from psydac.feec.global_projectors import Projector_H1 + from psydac.ddm.cart import DomainDecomposition from mpi4py import MPI @@ -348,13 +350,14 @@ def test_directional_derivative_operator_3d_par(domain, ncells, degree, periodic @pytest.mark.parametrize('degree', [2, 3, 4, 5]) @pytest.mark.parametrize('periodic', [True, False]) @pytest.mark.parametrize('seed', [1,3]) +@pytest.mark.parametrize('multiplicity', [1,2]) -def test_Derivative_1D(domain, ncells, degree, periodic, seed): +def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests np.random.seed(seed) breaks = np.linspace(*domain, num=ncells+1) - knots = make_knots(breaks, degree, periodic) + knots = make_knots(breaks, degree, periodic, multiplicity=multiplicity) # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') @@ -365,12 +368,11 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed): # L2 space (1-forms) V1 = V0.reduce_degree(axes=[0], basis='M') + u0 = FemField(V0) # Linear operator: 1D derivative grad = Derivative_1D(V0, V1) - + # Create random field in V0 - u0 = FemField(V0) - s, = V0.vector_space.starts e, = V0.vector_space.ends @@ -378,7 +380,7 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed): # Compute gradient (=derivative) of u0 u1 = grad(u0) - + #print(u1.coeffs) # Create evaluation grid, and check if ∂/∂x u0(x) == u1(x) xgrid = np.linspace(*N.domain, num=11) vals_grad_u0 = np.array([u0.gradient(x)[0] for x in xgrid]) @@ -395,8 +397,9 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed): @pytest.mark.parametrize('degree', [(3, 2), (4, 5)]) # 2 cases @pytest.mark.parametrize('periodic', [(True, False), (False, True)]) # 2 cases @pytest.mark.parametrize('seed', [1,3]) +@pytest.mark.parametrize('multiplicity', [[1,1],[1,2],[2,2]]) -def test_Gradient_2D(domain, ncells, degree, periodic, seed): +def test_Gradient_2D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests np.random.seed(seed) @@ -404,8 +407,8 @@ def test_Gradient_2D(domain, ncells, degree, periodic, seed): breaks = [np.linspace(*lims, num=n+1) for lims, n in zip(domain, ncells)] # H1 space (0-forms) - Nx, Ny = [SplineSpace(degree=d, grid=g, periodic=p, basis='B') \ - for d, g, p in zip(degree, breaks, periodic)] + Nx, Ny = [SplineSpace(degree=d, grid=g, periodic=p, basis='B', multiplicity=m) \ + for d, g, p, m in zip(degree, breaks, periodic, multiplicity)] domain_decomposition = DomainDecomposition(ncells, periodic) V0 = TensorFemSpace(domain_decomposition, Nx, Ny) @@ -456,11 +459,17 @@ def test_Gradient_2D(domain, ncells, degree, periodic, seed): (False, True, False), (False, False, True)]) @pytest.mark.parametrize('seed', [1,3]) +@pytest.mark.parametrize('multiplicity', [[1,1,1],[1,2,2],[2,2,2]]) -def test_Gradient_3D(domain, ncells, degree, periodic, seed): + +def test_Gradient_3D(domain, ncells, degree, periodic, seed, multiplicity): if any([ncells[d] <= degree[d] and periodic[d] for d in range(3)]): return + #change mulitplicity if higher than degree to avoid problems (case p Date: Fri, 20 Oct 2023 16:44:11 +0200 Subject: [PATCH 27/37] solved the problem, update tests and add tests in feec api --- psydac/api/fem.py | 4 +- psydac/api/grid.py | 6 +- psydac/api/tests/test_api_feec_1d.py | 8 +- psydac/api/tests/test_api_feec_2d.py | 24 +++--- psydac/api/tests/test_api_feec_3d.py | 2 +- psydac/core/bsplines.py | 42 ++++----- psydac/core/bsplines_kernels.py | 37 ++++---- psydac/core/tests/test_bsplines.py | 14 +-- psydac/core/tests/test_bsplines_pyccel.py | 67 +++++++-------- .../tests/test_differentiation_matrices.py | 59 +++++++------ psydac/feec/tests/test_global_projectors.py | 85 +++++++++++-------- psydac/fem/splines.py | 19 +++-- psydac/fem/tensor.py | 8 +- 13 files changed, 193 insertions(+), 182 deletions(-) diff --git a/psydac/api/fem.py b/psydac/api/fem.py index 20de7baaf..cf593d140 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -610,7 +610,7 @@ def construct_arguments(self, with_openmp=False): space = spaces[i].spaces[axis] points_i = points[i][axis] periodic = space.periodic - local_span = find_span(space.knots, space.degree, points_i[0, 0], periodic) + local_span = find_span(space.knots, space.degree, points_i[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points_i[0, 0], local_span, nderiv, space.basis) map_basis[i][axis] = map_basis[i][axis].copy() @@ -1214,7 +1214,7 @@ def construct_arguments(self, with_openmp=False): nderiv = self.max_nderiv space = space.spaces[axis] points = points[axis] - local_span = find_span(space.knots, space.degree, points[0, 0], space.periodic) + local_span = find_span(space.knots, space.degree, points[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) map_basis[axis] = map_basis[axis].copy() diff --git a/psydac/api/grid.py b/psydac/api/grid.py index 1771bc131..f04a1e01f 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -206,7 +206,7 @@ def __init__( self, V, nderiv , trial=False, grid=None): for i,Vi in enumerate(V): space = Vi.spaces[axis] points = grid.points[axis] - local_span = find_span(space.knots, space.degree, points[0, 0], space.periodic) + local_span = find_span(space.knots, space.degree, points[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points[0, 0], local_span, nderiv, space.basis) @@ -230,13 +230,11 @@ def space(self): # TODO have a parallel version of this function, as done for fem def create_collocation_basis( glob_points, space, nderiv=1 ): - #V.C 12/10/23 : Is this function really used somewhere? What is it suppose to do? T = space.knots # knots sequence p = space.degree # spline degree n = space.nbasis # total number of control points grid = space.breaks # breakpoints nc = space.ncells # number of cells in domain (nc=len(grid)-1) - perio = space.periodic # periodicity of the space #------------------------------------------- # GLOBAL GRID @@ -248,7 +246,7 @@ def create_collocation_basis( glob_points, space, nderiv=1 ): # glob_basis = np.zeros( (p+1,nderiv+1,nq) ) # TODO use this for local basis fct glob_basis = np.zeros( (n+p,nderiv+1,nq) ) # n+p for ghosts for iq,xq in enumerate(glob_points): - span = find_span( T, p, xq, perio ) + span = find_span( T, p, xq) glob_spans[iq] = span ders = basis_funs_all_ders( T, p, xq, span, nderiv ) diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index 93f69a45d..3949454db 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -415,8 +415,8 @@ def discrete_energies(self, e, b): errB = lambda x1: (push_1d_l2(B, x1, F) - B_ex(t, *F(x1)))**2 * np.sqrt(F.metric_det(x1)) error_l2_E = np.sqrt(derham_h.V1.integral(errE)) error_l2_B = np.sqrt(derham_h.V0.integral(errB)) - print('L2 norm of error on E(t,x) at final time: {:.2e}'.format(error_l2_E)) - print('L2 norm of error on B(t,x) at final time: {:.2e}'.format(error_l2_B)) + print('L2 norm of error on E(t,x) at final time: {:.12e}'.format(error_l2_E)) + print('L2 norm of error on B(t,x) at final time: {:.12e}'.format(error_l2_B)) if diagnostics_interval: @@ -518,8 +518,8 @@ def test_maxwell_1d_periodic_mult(): ) TOL = 1e-6 - ref = dict(error_E = 4.54259221e-04, - error_B = 3.71519383e-04) + ref = dict(error_E = 4.24689338e-04, + error_B = 4.03195792e-04) assert abs(namespace['error_E'] - ref['error_E']) / ref['error_E'] <= TOL assert abs(namespace['error_B'] - ref['error_B']) / ref['error_B'] <= TOL diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index c0c8f0c0d..2db063d02 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -730,9 +730,9 @@ def test_maxwell_2d_multiplicity(): ) TOL = 1e-6 - ref = dict(error_l2_Ex = 4.0298626142e-02, - error_l2_Ey = 4.0298626142e-02, - error_l2_Bz = 3.0892842064e-01) + ref = dict(error_l2_Ex = 4.35005708e-04, + error_l2_Ey = 4.35005708e-04, + error_l2_Bz = 3.76106860e-03) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL @@ -746,8 +746,7 @@ def test_maxwell_2d_periodic_multiplicity(): ncells = 30, degree = 3, periodic = True, - Cp = 0.5, #it's seems that with higher multiplicity we need - #to be more carefull with the CFL (with C=5. error is ) + Cp = 0.5, nsteps = 1, tend = None, splitting_order = 2, @@ -759,9 +758,9 @@ def test_maxwell_2d_periodic_multiplicity(): ) TOL = 1e-6 - ref = dict(error_l2_Ex = 8.2398553953e-02, - error_l2_Ey = 8.2398553953e-02, - error_l2_Bz = 2.6713526134e-01) + ref = dict(error_l2_Ex = 1.78557685e-04, + error_l2_Ey = 1.78557685e-04, + error_l2_Bz = 1.40582413e-04) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL @@ -788,9 +787,9 @@ def test_maxwell_2d_periodic_multiplicity_equal_deg(): ) TOL = 1e-6 - ref = dict(error_l2_Ex = 7.43229667e-02, - error_l2_Ey = 7.43229667e-02, - error_l2_Bz = 2.59863979e-01) + ref = dict(error_l2_Ex = 2.50585008e-02, + error_l2_Ey = 2.50585008e-02, + error_l2_Bz = 1.58438290e-02) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL assert abs(namespace['error_l2_Ey'] - ref['error_l2_Ey']) / ref['error_l2_Ey'] <= TOL @@ -913,9 +912,6 @@ def test_maxwell_2d_dirichlet_par(): # SCRIPT CAPABILITIES #============================================================================== if __name__ == '__main__': - - test_maxwell_2d_multiplicity() - exit() import argparse diff --git a/psydac/api/tests/test_api_feec_3d.py b/psydac/api/tests/test_api_feec_3d.py index 1bbf53318..ac7db93c0 100644 --- a/psydac/api/tests/test_api_feec_3d.py +++ b/psydac/api/tests/test_api_feec_3d.py @@ -395,7 +395,7 @@ class CollelaMapping3D(Mapping): T = dt*niter error = run_maxwell_3d_stencil(logical_domain, M, e_ex, b_ex, ncells, degree, periodic, dt, niter, mult=2) - assert abs(error - 0.9020674104318802) < 1e-9 + assert abs(error - 0.24749763720543216) < 1e-9 #============================================================================== # CLEAN UP SYMPY NAMESPACE diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 2fa85aa05..13fc22026 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -57,7 +57,7 @@ #============================================================================== -def find_span(knots, degree, x, periodic, multiplicity = 1): +def find_span(knots, degree, x): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -74,15 +74,7 @@ def find_span(knots, degree, x, periodic, multiplicity = 1): Polynomial degree of B-splines. x : float - Location of interest. - - periodic : bool - Indicating if the domain is periodic. - - multiplicity : int - Multplicity in the knot sequence. One assume that all the knots have - multiplicity multiplicity, and in the none periodic case that the - boundary knots have multiplicity p+1 (has created by make_knots) + Location of interest.. Returns ------- @@ -91,12 +83,10 @@ def find_span(knots, degree, x, periodic, multiplicity = 1): """ x = float(x) knots = np.ascontiguousarray(knots, dtype=float) - multiplicity = int(multiplicity) - periodic = bool(periodic) - return find_span_p(knots, degree, x, periodic, multiplicity = multiplicity) + return find_span_p(knots, degree, x) #============================================================================== -def find_spans(knots, degree, x, periodic, out=None): +def find_spans(knots, degree, x, out=None): """ Determine the knot span index at a set of locations x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -131,7 +121,7 @@ def find_spans(knots, degree, x, periodic, out=None): else: assert out.shape == x.shape and out.dtype == np.dtype('int') - find_spans_p(knots, degree, x, periodic, out) + find_spans_p(knots, degree, x, out) return out #============================================================================== @@ -336,6 +326,10 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, out : array, optional If provided, the result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. Returns ------- @@ -352,7 +346,7 @@ def collocation_matrix(knots, degree, periodic, normalization, xgrid, out=None, if out is None: nb = len(knots) - degree - 1 if periodic: - nb -= degree + nb -= degree + 1 - multiplicity out = np.zeros((xgrid.shape[0], nb), dtype=float) else: @@ -386,6 +380,10 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli xgrid : array_like Grid points. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. check_boundary : bool, default=True If true and ``periodic``, will check the boundaries of ``xgrid``. @@ -435,12 +433,12 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, multipli if out is None: if periodic: - out = np.zeros((len(xgrid), len(knots) - 2 * degree - 1), dtype=float) + out = np.zeros((len(xgrid), len(knots) - 2 * degree - 2 + multiplicity), dtype=float) else: out = np.zeros((len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1), dtype=float) else: if periodic: - assert out.shape == (len(xgrid), len(knots) - 2 * degree - 1) + assert out.shape == (len(xgrid), len(knots) - 2 * degree - 2 + multiplicity) else: assert out.shape == (len(xgrid) - 1, len(elevated_knots) - (degree + 1) - 1 - 1) assert out.dtype == np.dtype('float') @@ -501,6 +499,10 @@ def greville(knots, degree, periodic, out=None, multiplicity=1): out : array, optional If provided, the result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. Returns ------- @@ -510,7 +512,7 @@ def greville(knots, degree, periodic, out=None, multiplicity=1): """ knots = np.ascontiguousarray(knots, dtype=float) if out is None: - n = len(knots) - 2 * degree - 1 if periodic else len(knots) - degree - 1 + n = len(knots) - 2 * degree - 2 + multiplicity if periodic else len(knots) - degree - 1 out = np.zeros(n) multiplicity = int(multiplicity) greville_p(knots, degree, periodic, out, multiplicity) @@ -1011,7 +1013,7 @@ def alpha_function(i, k, t, n, p, knots): mat = np.zeros((n+1,n)) - left = find_span( knots, p, t, False ) + left = find_span( knots, p, t ) # ... j = 0 diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index d5c60be6a..b41d3305f 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -5,7 +5,7 @@ import numpy as np # ============================================================================= -def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multiplicity : int = 1): +def find_span_p(knots: 'float[:]', degree: int, x: float): """ Determine the knot span index at location x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -24,14 +24,6 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi x : float Location of interest. - periodic : bool - Indicating if the domain is periodic. - - multiplicity : int - Multplicity in the knot sequence. One assume that all the knots have - multiplicity multiplicity, and in the none periodic case that the - boundary knots have multiplicity p+1 (has created by make_knots) - Returns ------- span : int @@ -43,9 +35,6 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi Springer-Verlag Berlin Heidelberg GmbH, 1997. """ # Knot index at left/right boundary - #In the periodic case the p+1 knots is the last knot equal to the left - #boundary, while the first knots being the right boundary is multiplicity+p - #before the last knot of the sequence (see make knots) low = degree high = len(knots)-1-degree @@ -66,7 +55,7 @@ def find_span_p(knots: 'float[:]', degree: int, x: float, periodic : bool, multi # ============================================================================= -def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', periodic : bool, multiplicity : int = 1): +def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]'): """ Determine the knot span index at a set of locations x, given the B-Splines' knot sequence and polynomial degree. See Algorithm A2.1 in [1]. @@ -92,7 +81,7 @@ def find_spans_p(knots: 'float[:]', degree: int, x: 'float[:]', out: 'int[:]', p n = x.shape[0] for i in range(n): - out[i] = find_span_p(knots, degree, x[i], periodic, multiplicity=multiplicity) + out[i] = find_span_p(knots, degree, x[i]) # ============================================================================= @@ -422,11 +411,15 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali The result will be inserted into this array. It should be of the appropriate shape and dtype. + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. + """ # Number of basis functions (in periodic case remove degree repeated elements) nb = len(knots)-degree-1 if periodic: - nb -= degree + nb -= degree + 1 - multiplicity # Number of evaluation points nx = len(xgrid) @@ -434,7 +427,7 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali basis = np.zeros((nx, degree + 1)) spans = np.zeros(nx, dtype=int) - find_spans_p(knots, degree, xgrid, spans, periodic, multiplicity = multiplicity) + find_spans_p(knots, degree, xgrid, spans) basis_funs_array_p(knots, degree, xgrid, spans, basis) # Fill in non-zero matrix values @@ -500,6 +493,10 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma out : array The result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. Notes ----- @@ -509,7 +506,7 @@ def histopolation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, norma """ nb = len(knots) - degree - 1 if periodic: - nb -= degree + nb -= degree + 1 - multiplicity # Number of evaluation points nx = len(xgrid) @@ -718,10 +715,14 @@ def greville_p(knots: 'float[:]', degree: int, periodic: bool, out:'float[:]', m out : array The result will be inserted into this array. It should be of the appropriate shape and dtype. + + multiplicity : int + Multiplicity of the knots in the knot sequence, we assume that the same + multiplicity applies to each interior knot. """ T = knots p = degree - n = len(T)-2*p-1 if periodic else len(T)-p-1 + n = len(T)-2*p-2 + multiplicity if periodic else len(T)-p-1 # Compute greville abscissas as average of p consecutive knot values if p == multiplicity-1: diff --git a/psydac/core/tests/test_bsplines.py b/psydac/core/tests/test_bsplines.py index c539bbd9a..4b3f0b987 100644 --- a/psydac/core/tests/test_bsplines.py +++ b/psydac/core/tests/test_bsplines.py @@ -32,9 +32,9 @@ def test_find_span( lims, nc, p, eps=1e-12 ): knots = np.r_[ [grid[0]]*p, grid, [grid[-1]]*p ] for i,xi in enumerate( grid ): - assert find_span( knots, p, x=xi-eps, periodic=False ) == p + max( 0, i-1 ) - assert find_span( knots, p, x=xi, periodic=False ) == p + min( i, nc-1 ) - assert find_span( knots, p, x=xi+eps, periodic=False ) == p + min( i, nc-1 ) + assert find_span( knots, p, x=xi-eps ) == p + max( 0, i-1 ) + 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( 'lims', ([0,1], [-2,3]) ) @@ -48,7 +48,7 @@ def test_basis_funs( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) basis = basis_funs( knots, p, x, span ) assert len( basis ) == p+1 assert np.all( basis >= 0 ) @@ -66,7 +66,7 @@ def test_basis_funs_1st_der( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) ders = basis_funs_1st_der( knots, p, x, span ) assert len( ders ) == p+1 assert abs( sum( ders ) ) < tol @@ -86,7 +86,7 @@ def test_basis_funs_all_ders( lims, nc, p, tol=1e-14 ): xx = np.linspace( *lims, num=101 ) for x in xx: - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) ders = basis_funs_all_ders( knots, p, x, span, n ) # Test output array @@ -195,7 +195,7 @@ def test_cell_index(i_grid, expected): yy = np.zeros( (len(xx), nb) ) zz = np.zeros( (len(xx), nb) ) for i,x in enumerate( xx ): - span = find_span( knots, p, x, False ) + span = find_span( knots, p, x ) yy[i,span-p:span+1] = basis_funs ( knots, p, x, span ) zz[i,span-p:span+1] = basis_funs_1st_der( knots, p, x, span ) diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index aeb9f68cd..f243f39ae 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -21,13 +21,11 @@ # "True" Functions ############################################################################### -def find_span_true( knots, degree, x, periodic, multiplicity=1 ): +def find_span_true( knots, degree, x): # Knot index at left/right boundary low = degree - if periodic : - high = len(knots)-multiplicity-degree - else : - high = len(knots)-1-degree + high = len(knots)-1-degree + # Check if point is exactly on left/right boundary, or outside domain if x <= knots[low ]: return low if x >= knots[high]: return high-1 @@ -211,7 +209,7 @@ def collocation_matrix_true(knots, degree, periodic, normalization, xgrid): # Fill in non-zero matrix values for i,x in enumerate( xgrid ): - span = find_span_true( knots, degree, x, periodic ) + span = find_span( knots, degree, x ) basis = basis_funs_true( knots, degree, x, span ) mat[i,js(span)] = normalize(basis, span) @@ -377,27 +375,29 @@ def make_knots_true( breaks, degree, periodic, multiplicity=1 ): if periodic: assert len(breaks) > degree - - p = degree - if periodic : - T = np.zeros(multiplicity * len(breaks[1:]) + 1 + 2 * degree) - else : - T = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) - + + T = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) + + ncells = len(breaks) - 1 + + for i in range(0, ncells+1): + T[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] + + len_out = len(T) + if periodic: period = breaks[-1]-breaks[0] - T[p+1:-p] = np.repeat(breaks[1:], multiplicity) - T[0:p+1] = [xi-period for xi in T[-2*p-1:-p ]] - T[-p:] = [xi+period for xi in T[ p+1:2*p+1]] + T[: degree + 1 - multiplicity] = T[len_out - 2 * (degree + 1 )+ multiplicity: len_out - (degree + 1)] - period + T[len_out - (degree + 1 - multiplicity) :] = T[degree + 1:2*(degree + 1)- multiplicity] + period + + + # else: - T[p+1:-p-1] = np.repeat(breaks[1:-1], multiplicity) - T[p] = breaks[ 0] - T[-p-1] = breaks[-1] - T[0:p+1] = breaks[ 0] - T[-p-1:] = breaks[-1] + T[0:degree + 1 - multiplicity] = breaks[0] + T[len_out - degree - 1 + multiplicity:] = breaks[-1] return T @@ -406,10 +406,10 @@ def elevate_knots_true(knots, degree, periodic, multiplicity=1, tol=1e-15): knots = np.array(knots) if periodic: - [T, p] = knots, degree - period = T[-1-p] - T[p] - left = [T[-1-p-(p+1)] - period] - right = [T[ p+(p+1)] + period] + T, p = knots, degree + period = T[len(knots) -1 - p] - T[p] + left = [T[len(knots) -2 - 2 * p + multiplicity-1] - period] + right = [T[2 * p + 2 - multiplicity] + period] else: left = [knots[0],*knots[:degree+1]] right = [knots[-1],*knots[-degree-1:]] @@ -461,7 +461,7 @@ def basis_ders_on_quad_grid_true(knots, degree, quad_grid, nders, normalization) for ie in range(ne): xx = quad_grid[ie, :] for iq, xq in enumerate(xx): - span = find_span_true(knots, degree, xq, False) + span = find_span_true(knots, degree, xq) ders = basis_funs_all_ders_true(knots, degree, xq, span, nders) if normalization == 'M': ders *= scaling[None, span-degree:span+1] @@ -499,11 +499,10 @@ def basis_integrals_true(knots, degree): (np.array([0.0, 0.0, 0.0, 1.0, 1.0, 1.0]), 2), (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) -@pytest.mark.parametrize('periodic', [True, False]) -@pytest.mark.parametrize('multiplicity', [1, 2, 3]) -def test_find_span(knots, degree, x, periodic, multiplicity): - expected = find_span_true(knots, degree, x, periodic, multiplicity=multiplicity) - out = find_span(knots, degree, x, periodic, multiplicity=multiplicity) + +def test_find_span(knots, degree, x): + expected = find_span_true(knots, degree, x) + out = find_span(knots, degree, x) assert np.allclose(expected, out, atol=ATOL, rtol=RTOL) @@ -521,7 +520,7 @@ def test_find_span(knots, degree, x, periodic, multiplicity): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs(knots, degree, x): - span = find_span(knots, degree, x, False) + span = find_span(knots, degree, x) expected = basis_funs_true(knots, degree, x, span) out = basis_funs(knots, degree, x, span) @@ -541,7 +540,7 @@ def test_basis_funs(knots, degree, x): (np.array([0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]), 3)]) @pytest.mark.parametrize('x', (np.random.random(), np.random.random(), np.random.random())) def test_basis_funs_1st_der(knots, degree, x): - span = find_span(knots, degree, x, False) + span = find_span(knots, degree, x) expected = basis_funs_1st_der_true(knots, degree, x, span) out = basis_funs_1st_der(knots, degree, x, span) @@ -563,7 +562,7 @@ def test_basis_funs_1st_der(knots, degree, x): @pytest.mark.parametrize('n', (2, 3, 4, 5)) @pytest.mark.parametrize('normalization', ('B', 'M')) def test_basis_funs_all_ders(knots, degree, x, n, normalization): - span = find_span(knots, degree, x, False) + span = find_span(knots, degree, x) expected = basis_funs_all_ders_true(knots, degree, x, span, n, normalization) out = basis_funs_all_ders(knots, degree, x, span, n, normalization) diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 3feccd55d..35e017421 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -461,21 +461,18 @@ def test_Gradient_2D(domain, ncells, degree, periodic, seed, multiplicity): @pytest.mark.parametrize('seed', [1,3]) @pytest.mark.parametrize('multiplicity', [[1,1,1],[1,2,2],[2,2,2]]) - def test_Gradient_3D(domain, ncells, degree, periodic, seed, multiplicity): if any([ncells[d] <= degree[d] and periodic[d] for d in range(3)]): return - #change mulitplicity if higher than degree to avoid problems (case p1e-15)[0] - if len(indices)>0: - multiplicity = np.diff(indices).max(initial=1) - else: - multiplicity = max(1,len(knots[degree+1:-degree-1])) + if multiplicity is None: #V.C 20/10/23 : why computing the multiplicity again? + if len(indices)>0: + multiplicity = np.diff(indices).max(initial=1) + else: + multiplicity = max(1,len(knots[degree+1:-degree-1])) + if parent_multiplicity is None: parent_multiplicity = multiplicity @@ -107,7 +109,7 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul # Number of basis function in space (= cardinality) if periodic: - nbasis = len(knots) - 2*degree - 1 + nbasis = len(knots) - 2*degree - 2 + multiplicity else: defect = 0 if dirichlet[0]: defect += 1 @@ -134,7 +136,6 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul self._greville = greville(knots, degree, periodic, multiplicity = multiplicity) self._ext_greville = greville(elevate_knots(knots, degree, periodic, multiplicity=multiplicity), degree+1, periodic, multiplicity = multiplicity) self._scaling_array = scaling_array - self._parent_multiplicity = parent_multiplicity self._histopolation_grid = unroll_edges(self.domain, self.ext_greville) @@ -181,7 +182,7 @@ def init_interpolation( self, dtype=float ): xgrid = self.greville, multiplicity = self.multiplicity ) - + if self.periodic: # Convert to CSC format and compute sparse LU decomposition self._interpolator = SparseSolver( csc_matrix( imat ) ) @@ -290,7 +291,7 @@ def eval_field(self, field, *eta , weights=None): eta = eta[0] - span = find_span( self.knots, self.degree, eta, self.periodic) + span = find_span( self.knots, self.degree, eta) basis_array = basis_funs( self.knots, self.degree, eta, span) index = slice(span-self.degree, span + 1) diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 6a9b0b86b..4ea6a9b74 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -193,9 +193,7 @@ def eval_field( self, field, *eta, weights=None): for (x, xlim, space) in zip( eta, self.eta_lims, self.spaces ): knots = space.knots degree = space.degree - multiplicity = space.multiplicity - periodic = space.periodic - span = find_span( knots, degree, x, periodic, multiplicity = multiplicity) + span = find_span( knots, degree, x ) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# @@ -622,9 +620,7 @@ def eval_field_gradient( self, field, *eta , weights=None): knots = space.knots degree = space.degree - multiplicity = space.multiplicity - periodic = space.periodic - span = find_span( knots, degree, x, periodic, multiplicity = multiplicity) + span = find_span( knots, degree, x ) #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# From 4c14de1017530a787965ff76be759e58ec5090c7 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Mon, 23 Oct 2023 11:02:19 +0200 Subject: [PATCH 28/37] clean the PR --- psydac/api/fem.py | 1 - psydac/api/grid.py | 14 +++++++------- psydac/api/tests/test_api_feec_1d.py | 8 ++++---- psydac/api/tests/test_assembly.py | 2 +- psydac/core/bsplines.py | 2 +- psydac/core/bsplines_kernels.py | 6 +++--- psydac/core/tests/test_bsplines_pyccel.py | 4 ++-- psydac/feec/tests/test_differentiation_matrices.py | 3 +-- psydac/fem/splines.py | 2 +- psydac/fem/tensor.py | 2 ++ 10 files changed, 22 insertions(+), 22 deletions(-) diff --git a/psydac/api/fem.py b/psydac/api/fem.py index cf593d140..c96bc9b56 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -609,7 +609,6 @@ def construct_arguments(self, with_openmp=False): if axis is None:continue space = spaces[i].spaces[axis] points_i = points[i][axis] - periodic = space.periodic local_span = find_span(space.knots, space.degree, points_i[0, 0]) boundary_basis = basis_funs_all_ders(space.knots, space.degree, points_i[0, 0], local_span, nderiv, space.basis) diff --git a/psydac/api/grid.py b/psydac/api/grid.py index f04a1e01f..6e8b3f118 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -229,12 +229,12 @@ def space(self): #============================================================================== # TODO have a parallel version of this function, as done for fem def create_collocation_basis( glob_points, space, nderiv=1 ): - - T = space.knots # knots sequence - p = space.degree # spline degree - n = space.nbasis # total number of control points - grid = space.breaks # breakpoints - nc = space.ncells # number of cells in domain (nc=len(grid)-1) + + T = space.knots # knots sequence + p = space.degree # spline degree + n = space.nbasis # total number of control points + grid = space.breaks # breakpoints + nc = space.ncells # number of cells in domain (nc=len(grid)-1) #------------------------------------------- # GLOBAL GRID @@ -246,7 +246,7 @@ def create_collocation_basis( glob_points, space, nderiv=1 ): # glob_basis = np.zeros( (p+1,nderiv+1,nq) ) # TODO use this for local basis fct glob_basis = np.zeros( (n+p,nderiv+1,nq) ) # n+p for ghosts for iq,xq in enumerate(glob_points): - span = find_span( T, p, xq) + span = find_span( T, p, xq ) glob_spans[iq] = span ders = basis_funs_all_ders( T, p, xq, span, nderiv ) diff --git a/psydac/api/tests/test_api_feec_1d.py b/psydac/api/tests/test_api_feec_1d.py index 3949454db..260bb0381 100644 --- a/psydac/api/tests/test_api_feec_1d.py +++ b/psydac/api/tests/test_api_feec_1d.py @@ -406,8 +406,8 @@ def discrete_energies(self, e, b): error_E = max(abs(E_ex(t, x) - E_values)) error_B = max(abs(B_ex(t, x) - B_values)) print() - print('Max-norm of error on E(t,x) at final time: {:.8e}'.format(error_E)) - print('Max-norm of error on B(t,x) at final time: {:.8e}'.format(error_B)) + print('Max-norm of error on E(t,x) at final time: {:.2e}'.format(error_E)) + print('Max-norm of error on B(t,x) at final time: {:.2e}'.format(error_B)) # compute L2 error as well F = mapping.get_callable_mapping() @@ -415,8 +415,8 @@ def discrete_energies(self, e, b): errB = lambda x1: (push_1d_l2(B, x1, F) - B_ex(t, *F(x1)))**2 * np.sqrt(F.metric_det(x1)) error_l2_E = np.sqrt(derham_h.V1.integral(errE)) error_l2_B = np.sqrt(derham_h.V0.integral(errB)) - print('L2 norm of error on E(t,x) at final time: {:.12e}'.format(error_l2_E)) - print('L2 norm of error on B(t,x) at final time: {:.12e}'.format(error_l2_B)) + print('L2 norm of error on E(t,x) at final time: {:.2e}'.format(error_l2_E)) + print('L2 norm of error on B(t,x) at final time: {:.2e}'.format(error_l2_B)) if diagnostics_interval: diff --git a/psydac/api/tests/test_assembly.py b/psydac/api/tests/test_assembly.py index b598eb55d..7fdbdfce1 100644 --- a/psydac/api/tests/test_assembly.py +++ b/psydac/api/tests/test_assembly.py @@ -117,7 +117,7 @@ def test_field_and_constant_mult(backend, dtype): ncells = (5, 5) degree = (3, 3) - multiplicity = (2,2) + multiplicity = [2,2] domain_h = discretize(domain, ncells=ncells) Vh = discretize(V, domain_h, degree=degree, multiplicity = multiplicity) ah = discretize(a, domain_h, [Vh, Vh], **kwargs) diff --git a/psydac/core/bsplines.py b/psydac/core/bsplines.py index 13fc22026..3567b1f1a 100644 --- a/psydac/core/bsplines.py +++ b/psydac/core/bsplines.py @@ -74,7 +74,7 @@ def find_span(knots, degree, x): Polynomial degree of B-splines. x : float - Location of interest.. + Location of interest. Returns ------- diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index b41d3305f..f4059e8c0 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -23,7 +23,7 @@ def find_span_p(knots: 'float[:]', degree: int, x: float): x : float Location of interest. - + Returns ------- span : int @@ -426,11 +426,11 @@ def collocation_matrix_p(knots: 'float[:]', degree: int, periodic: bool, normali basis = np.zeros((nx, degree + 1)) spans = np.zeros(nx, dtype=int) - find_spans_p(knots, degree, xgrid, spans) basis_funs_array_p(knots, degree, xgrid, spans, basis) # Fill in non-zero matrix values + # Rescaling of B-splines, to get M-splines if needed if not normalization: if periodic: @@ -838,7 +838,7 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] period = breaks[-1]-breaks[0] out[: degree + 1 - multiplicity] = out[len_out - 2 * (degree + 1 )+ multiplicity: len_out - (degree + 1)] - period - out[len(out) - (degree + 1 - multiplicity) :] = out[degree + 1:2*(degree + 1)- multiplicity] + period + out[len_out - (degree + 1 - multiplicity) :] = out[degree + 1:2*(degree + 1)- multiplicity] + period # diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index f243f39ae..a1193b94d 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -21,7 +21,7 @@ # "True" Functions ############################################################################### -def find_span_true( knots, degree, x): +def find_span_true( knots, degree, x ): # Knot index at left/right boundary low = degree high = len(knots)-1-degree @@ -209,7 +209,7 @@ def collocation_matrix_true(knots, degree, periodic, normalization, xgrid): # Fill in non-zero matrix values for i,x in enumerate( xgrid ): - span = find_span( knots, degree, x ) + span = find_span_true( knots, degree, x ) basis = basis_funs_true( knots, degree, x, span ) mat[i,js(span)] = normalize(basis, span) diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 35e017421..6d56cc6c9 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -371,7 +371,7 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): u0 = FemField(V0) # Linear operator: 1D derivative grad = Derivative_1D(V0, V1) - + # Create random field in V0 s, = V0.vector_space.starts e, = V0.vector_space.ends @@ -380,7 +380,6 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): # Compute gradient (=derivative) of u0 u1 = grad(u0) - #print(u1.coeffs) # Create evaluation grid, and check if ∂/∂x u0(x) == u1(x) xgrid = np.linspace(*N.domain, num=11) vals_grad_u0 = np.array([u0.gradient(x)[0] for x in xgrid]) diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 628d7f104..bfd7d9c6a 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -90,7 +90,7 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul if grid is None: grid = breakpoints(knots, degree) - + indices = np.where(np.diff(knots[degree+1:-degree-1])>1e-15)[0] diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 4ea6a9b74..e2e363e2d 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -191,9 +191,11 @@ 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 ): + knots = space.knots degree = space.degree span = find_span( knots, degree, x ) + #-------------------------------------------------# # Fix span for boundaries between subdomains # #-------------------------------------------------# From ab7c6c06c5a42c61a8d2a0af032182372041a8f0 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Mon, 23 Oct 2023 11:10:37 +0200 Subject: [PATCH 29/37] remove test in test_assembly --- psydac/api/tests/test_assembly.py | 56 ------------------------------- 1 file changed, 56 deletions(-) diff --git a/psydac/api/tests/test_assembly.py b/psydac/api/tests/test_assembly.py index 7fdbdfce1..4c52eeca7 100644 --- a/psydac/api/tests/test_assembly.py +++ b/psydac/api/tests/test_assembly.py @@ -85,60 +85,6 @@ def test_field_and_constant(backend, dtype): assert abs(b.toarray().sum() - res) < 1e-12 print("PASSED") -#============================================================================== -def test_field_and_constant_mult(backend, dtype): - - # If 'backend' is specified, accelerate Python code by passing **kwargs - # to discretization of bilinear forms, linear forms and functionals. - kwargs = {'backend': PSYDAC_BACKENDS[backend]} if backend else {} - - domain = Square() - V = ScalarFunctionSpace('V', domain) - - # TODO: remove codomain_type when It is implemented in sympde - u = element_of(V, name='u') - v = element_of(V, name='v') - f = element_of(V, name='f') - - if dtype == 'complex': - c = Constant(name='c', complex=True) - V.codomain_type = dtype - g = I * c * f**2 - res = 1.j - cst=complex(1.0) - else: - c = Constant(name='c', real=True) - g = c * f**2 - res = 1 - cst=1.0 - - a = BilinearForm((u, v), integral(domain, g * u * v)) - l = LinearForm(v, integral(domain, g * v)) - - ncells = (5, 5) - degree = (3, 3) - multiplicity = [2,2] - domain_h = discretize(domain, ncells=ncells) - Vh = discretize(V, domain_h, degree=degree, multiplicity = multiplicity) - ah = discretize(a, domain_h, [Vh, Vh], **kwargs) - lh = discretize(l, domain_h, Vh , **kwargs) - - fh = FemField(Vh) - fh.coeffs[:] = 1 - - # Assembly call should not crash if correct arguments are used - A = ah.assemble(c=cst, f=fh) - b = lh.assemble(f=fh, c=cst) - - # Test matrix A - x = fh.coeffs - #TODO change res into np.conj(res) when the conjugate is applied in the dot product in sympde - assert abs(x.dot(A.dot(x)) - res) < 1e-12 - - # Test vector b - assert abs(b.toarray().sum() - res) < 1e-12 - print("PASSED") - #============================================================================== def test_bilinearForm_complex(backend): @@ -586,8 +532,6 @@ def test_assembly_no_synchr_args(backend): #============================================================================== if __name__ == '__main__': - test_field_and_constant_mult('pyccel-gcc','real') - exit() test_field_and_constant(None,'real') test_multiple_fields(None) test_math_imports(None) From f10d6330643f557e00342a18ec56c4d766d7dd27 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Mon, 23 Oct 2023 11:13:05 +0200 Subject: [PATCH 30/37] remove some stuff --- psydac/api/tests/test_assembly.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/api/tests/test_assembly.py b/psydac/api/tests/test_assembly.py index 4c52eeca7..0dfae99f5 100644 --- a/psydac/api/tests/test_assembly.py +++ b/psydac/api/tests/test_assembly.py @@ -532,7 +532,7 @@ def test_assembly_no_synchr_args(backend): #============================================================================== if __name__ == '__main__': - test_field_and_constant(None,'real') + test_field_and_constant(None) test_multiple_fields(None) test_math_imports(None) test_non_symmetric_BilinearForm(None) From 083009a58529083cdb1726279878ac80b030e589 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Mon, 23 Oct 2023 13:30:48 +0200 Subject: [PATCH 31/37] fix computation of multiplicity --- psydac/api/tests/test_api_feec_2d.py | 6 +++--- psydac/feec/tests/test_differentiation_matrices.py | 5 +++-- psydac/fem/splines.py | 12 +++++------- 3 files changed, 11 insertions(+), 12 deletions(-) diff --git a/psydac/api/tests/test_api_feec_2d.py b/psydac/api/tests/test_api_feec_2d.py index 2db063d02..5c27b678e 100644 --- a/psydac/api/tests/test_api_feec_2d.py +++ b/psydac/api/tests/test_api_feec_2d.py @@ -729,9 +729,9 @@ def test_maxwell_2d_multiplicity(): mult = 2 ) - TOL = 1e-6 - ref = dict(error_l2_Ex = 4.35005708e-04, - error_l2_Ey = 4.35005708e-04, + TOL = 1e-5 + ref = dict(error_l2_Ex = 4.350041934920621e-04, + error_l2_Ey = 4.350041934920621e-04, error_l2_Bz = 3.76106860e-03) assert abs(namespace['error_l2_Ex'] - ref['error_l2_Ex']) / ref['error_l2_Ex'] <= TOL diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 6d56cc6c9..252b3e192 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -911,10 +911,11 @@ def eval_div(fx, fy, fz, *eta): #============================================================================== if __name__ == '__main__': - - test_Derivative_1D(domain=[0, 1], ncells=12, degree=3, periodic=False, seed=1, multiplicity=1) + + test_Derivative_1D(domain=[0, 1], ncells=3, degree=3, periodic=False, seed=1, multiplicity=1) test_Derivative_1D(domain=[0, 1], ncells=12, degree=3, periodic=True, seed=1, multiplicity=1) + test_Gradient_2D( domain = ([0, 1], [0, 1]), ncells = (10, 15), diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index bfd7d9c6a..3612b51d0 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -91,15 +91,13 @@ def __init__( self, degree, knots=None, grid=None, multiplicity=None, parent_mul if grid is None: grid = breakpoints(knots, degree) - indices = np.where(np.diff(knots[degree+1:-degree-1])>1e-15)[0] + indices = np.where(np.diff(knots[degree:len(knots)-degree])>1e-15)[0] + if len(indices)>0: + multiplicity = np.diff(indices).max(initial=1) + else: + multiplicity = max(1,len(knots[degree+1:-degree-1])) - if multiplicity is None: #V.C 20/10/23 : why computing the multiplicity again? - if len(indices)>0: - multiplicity = np.diff(indices).max(initial=1) - else: - multiplicity = max(1,len(knots[degree+1:-degree-1])) - if parent_multiplicity is None: parent_multiplicity = multiplicity From 49df6ee902e71ad2b3d65aee865fbd914d17d9b4 Mon Sep 17 00:00:00 2001 From: vcarlier Date: Mon, 23 Oct 2023 14:03:50 +0200 Subject: [PATCH 32/37] update tests --- psydac/feec/tests/test_differentiation_matrices.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 252b3e192..dfd1b3360 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -911,7 +911,7 @@ def eval_div(fx, fy, fz, *eta): #============================================================================== if __name__ == '__main__': - + test_Derivative_1D(domain=[0, 1], ncells=3, degree=3, periodic=False, seed=1, multiplicity=1) test_Derivative_1D(domain=[0, 1], ncells=12, degree=3, periodic=True, seed=1, multiplicity=1) From 08ede9c7ba41cd22a52c4aa5cbdb68138c93403d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 23 Oct 2023 18:37:04 +0200 Subject: [PATCH 33/37] Apply suggestions from code review: cosmetic changes --- psydac/feec/tests/test_differentiation_matrices.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index dfd1b3360..a2a4f4cb4 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -396,7 +396,7 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): @pytest.mark.parametrize('degree', [(3, 2), (4, 5)]) # 2 cases @pytest.mark.parametrize('periodic', [(True, False), (False, True)]) # 2 cases @pytest.mark.parametrize('seed', [1,3]) -@pytest.mark.parametrize('multiplicity', [[1,1],[1,2],[2,2]]) +@pytest.mark.parametrize('multiplicity', [(1, 1), (1, 2), (2, 2)]) def test_Gradient_2D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests @@ -458,7 +458,7 @@ def test_Gradient_2D(domain, ncells, degree, periodic, seed, multiplicity): (False, True, False), (False, False, True)]) @pytest.mark.parametrize('seed', [1,3]) -@pytest.mark.parametrize('multiplicity', [[1,1,1],[1,2,2],[2,2,2]]) +@pytest.mark.parametrize('multiplicity', [(1, 1, 1), (1, 2, 2), (2, 2, 2)]) def test_Gradient_3D(domain, ncells, degree, periodic, seed, multiplicity): if any([ncells[d] <= degree[d] and periodic[d] for d in range(3)]): @@ -524,7 +524,7 @@ def test_Gradient_3D(domain, ncells, degree, periodic, seed, multiplicity): @pytest.mark.parametrize('degree', [(3, 2), (4, 5)]) # 2 cases @pytest.mark.parametrize('periodic', [(True, False), (False, True)]) # 2 cases @pytest.mark.parametrize('seed', [1,3]) -@pytest.mark.parametrize('multiplicity', [[1,1],[1,2],[2,2]]) +@pytest.mark.parametrize('multiplicity', [(1, 1), (1, 2), (2, 2)]) def test_ScalarCurl_2D(domain, ncells, degree, periodic, seed, multiplicity): @@ -600,7 +600,7 @@ def eval_curl(fx, fy, *eta): @pytest.mark.parametrize('degree', [(3, 2), (4, 5)]) # 2 cases @pytest.mark.parametrize('periodic', [(True, False), (False, True)]) # 2 cases @pytest.mark.parametrize('seed', [1,3]) -@pytest.mark.parametrize('multiplicity', [[1,1],[1,2],[2,2]]) +@pytest.mark.parametrize('multiplicity', [(1, 1), (1, 2), (2, 2)]) def test_VectorCurl_2D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests @@ -666,7 +666,7 @@ def eval_curl(f, *eta): (False, True, False), (False, False, True)]) @pytest.mark.parametrize('seed', [1,3]) -@pytest.mark.parametrize('multiplicity', [[1,1,1],[1,2,2],[2,2,2]]) +@pytest.mark.parametrize('multiplicity', [(1, 1, 1), (1, 2, 2), (2, 2, 2)]) def test_Curl_3D(domain, ncells, degree, periodic, seed, multiplicity): if any([ncells[d] <= degree[d] and periodic[d] for d in range(3)]): @@ -758,7 +758,7 @@ def eval_curl(fx, fy, fz, *eta): @pytest.mark.parametrize('degree', [(3, 2), (4, 5)]) # 2 cases @pytest.mark.parametrize('periodic', [(True, False), (False, True)]) # 2 cases @pytest.mark.parametrize('seed', [1,3]) -@pytest.mark.parametrize('multiplicity', [[1,1],[1,2],[2,2]]) +@pytest.mark.parametrize('multiplicity', [(1, 1), (1, 2), (2, 2)]) def test_Divergence_2D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests @@ -834,7 +834,7 @@ def eval_div(fx, fy, *eta): (False, True, False), (False, False, True)]) @pytest.mark.parametrize('seed', [1,3]) -@pytest.mark.parametrize('multiplicity', [[1,1,1],[1,2,2],[2,2,2]]) +@pytest.mark.parametrize('multiplicity', [(1, 1, 1), (1, 2, 2), (2, 2, 2)]) def test_Divergence_3D(domain, ncells, degree, periodic, seed, multiplicity): # determinize tests From 6cc6311b260ecccd815054ddb99000936c31f0ef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 23 Oct 2023 18:40:16 +0200 Subject: [PATCH 34/37] Apply suggestions from code review: cosmetic changes --- psydac/feec/tests/test_global_projectors.py | 32 +++++++++++---------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index 628941d40..520a9bed9 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -17,7 +17,7 @@ @pytest.mark.parametrize('ncells', [500]) @pytest.mark.parametrize('degree', [1, 2, 3, 4, 5, 6, 7]) @pytest.mark.parametrize('periodic', [False, True]) -@pytest.mark.parametrize('multiplicity', [1,2]) +@pytest.mark.parametrize('multiplicity', [1, 2]) def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity): @@ -55,9 +55,9 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity): @pytest.mark.parametrize('domain', [(0, 2*np.pi)]) @pytest.mark.parametrize('ncells', [100, 200, 300]) @pytest.mark.parametrize('degree', [2]) -@pytest.mark.parametrize('periodic', [True,False]) +@pytest.mark.parametrize('periodic', [True, False]) @pytest.mark.parametrize('nquads', [100, 120, 140, 160]) -@pytest.mark.parametrize('multiplicity', [1,2]) +@pytest.mark.parametrize('multiplicity', [1, 2]) def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity): @@ -69,8 +69,8 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity) domain_decomposition = DomainDecomposition([ncells], [periodic]) # H1 space (0-forms) - #change mulitplicity if higher than degree to avoid problems (case p Date: Tue, 24 Oct 2023 09:18:32 +0200 Subject: [PATCH 35/37] apply chnages required by Yaman --- psydac/core/tests/test_bsplines_pyccel.py | 5 ----- psydac/feec/global_projectors.py | 2 +- psydac/feec/tests/test_differentiation_matrices.py | 4 +++- psydac/feec/tests/test_global_projectors.py | 4 ++-- psydac/fem/splines.py | 2 +- 5 files changed, 7 insertions(+), 10 deletions(-) diff --git a/psydac/core/tests/test_bsplines_pyccel.py b/psydac/core/tests/test_bsplines_pyccel.py index a1193b94d..5d9cdaad2 100644 --- a/psydac/core/tests/test_bsplines_pyccel.py +++ b/psydac/core/tests/test_bsplines_pyccel.py @@ -377,9 +377,6 @@ def make_knots_true( breaks, degree, periodic, multiplicity=1 ): assert len(breaks) > degree T = np.zeros(multiplicity * len(breaks[1:-1]) + 2 + 2 * degree) - - - ncells = len(breaks) - 1 for i in range(0, ncells+1): @@ -392,9 +389,7 @@ def make_knots_true( breaks, degree, periodic, multiplicity=1 ): T[: degree + 1 - multiplicity] = T[len_out - 2 * (degree + 1 )+ multiplicity: len_out - (degree + 1)] - period T[len_out - (degree + 1 - multiplicity) :] = T[degree + 1:2*(degree + 1)- multiplicity] + period - - # else: T[0:degree + 1 - multiplicity] = breaks[0] T[len_out - degree - 1 + multiplicity:] = breaks[-1] diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index c651dd610..32ed88ef9 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -170,7 +170,7 @@ def __init__(self, space, nquads = None): solverblocks += [KroneckerLinearSolver(tensorspaces[i].vector_space, solvercells)] - dataslice = tuple(slice(p*m, -p*m) for p,m in zip(tensorspaces[i].vector_space.pads,tensorspaces[i].vector_space.shifts)) + dataslice = tuple(slice(p*m, -p*m) for p, m in zip(tensorspaces[i].vector_space.pads,tensorspaces[i].vector_space.shifts)) dofs[i] = rhsblocks[i]._data[dataslice] # finish arguments and create a lambda diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index dfd1b3360..0173147ba 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -369,6 +369,7 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): V1 = V0.reduce_degree(axes=[0], basis='M') u0 = FemField(V0) + # Linear operator: 1D derivative grad = Derivative_1D(V0, V1) @@ -380,6 +381,7 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): # Compute gradient (=derivative) of u0 u1 = grad(u0) + # Create evaluation grid, and check if ∂/∂x u0(x) == u1(x) xgrid = np.linspace(*N.domain, num=11) vals_grad_u0 = np.array([u0.gradient(x)[0] for x in xgrid]) @@ -470,7 +472,7 @@ def test_Gradient_3D(domain, ncells, degree, periodic, seed, multiplicity): # Compute breakpoints along each direction breaks = [np.linspace(*lims, num=n+1) for lims, n in zip(domain, ncells)] - #change mulitplicity if higher than degree to avoid problems (case p Date: Tue, 24 Oct 2023 09:27:19 +0200 Subject: [PATCH 36/37] cosmetic changes --- psydac/core/bsplines_kernels.py | 2 +- psydac/feec/tests/test_differentiation_matrices.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index f4059e8c0..725138404 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -833,7 +833,7 @@ def make_knots_p(breaks: 'float[:]', degree: int, periodic: bool, out: 'float[:] out[degree + 1 + (i-1) * multiplicity :degree + 1 + i * multiplicity ] = breaks[i] len_out = len(out) - + if periodic: period = breaks[-1]-breaks[0] diff --git a/psydac/feec/tests/test_differentiation_matrices.py b/psydac/feec/tests/test_differentiation_matrices.py index 7dc45c313..10fdf699a 100644 --- a/psydac/feec/tests/test_differentiation_matrices.py +++ b/psydac/feec/tests/test_differentiation_matrices.py @@ -368,6 +368,7 @@ def test_Derivative_1D(domain, ncells, degree, periodic, seed, multiplicity): # L2 space (1-forms) V1 = V0.reduce_degree(axes=[0], basis='M') + # Create random field in V0 u0 = FemField(V0) # Linear operator: 1D derivative From 464f5209b86315f22191bd595ee26eb41a6ede7d Mon Sep 17 00:00:00 2001 From: vcarlier Date: Wed, 25 Oct 2023 09:29:42 +0200 Subject: [PATCH 37/37] apply change requested by Martin --- psydac/core/bsplines_kernels.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index 725138404..a2d9baf51 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -34,8 +34,9 @@ def find_span_p(knots: 'float[:]', degree: int, x: float): .. [1] L. Piegl and W. Tiller. The NURBS Book, 2nd ed., Springer-Verlag Berlin Heidelberg GmbH, 1997. """ - # Knot index at left/right boundary + # last knot on the left boundary low = degree + # first knot on the right boundary high = len(knots)-1-degree # Check if point is exactly on left/right boundary, or outside domain