solved the problem, update tests and add tests in feec api

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

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 )

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

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
 

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

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

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,19 +411,23 @@ 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)
 
     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: