De Rham sequences w/ knot multiplicity > 1 (#339)

Enable the use of the de Rham sequence with arbitrary knots
multiplicity. With this we mean that the interior knots (all knots in
the periodic case) can have a constant multiplicity higher than 1.

List of modifications
---------------------
* Changes in `core/bsplines.py` and `core/bsplines_kernels.py` : 
  - Modify the `make_knots` function to create an appropriate knot
sequence in all cases. In the non periodic case the interior knots are
repeated `multiplicity` times and the boundary knots `p+1` times. In the
periodic case, interior knots are repeated `multiplicity` times and we
add `p+1` knots at each side by periodicity. In all cases the number of
knots is `(n_cells-1) * multiplicity + 2 * (p+1)`.
  - Modify functions `collocation_matrix` and `histopolation_matrix` to
take in account that now in the periodic case the number of basis
function is `len(knots) - 2 * (p+1) + multiplicity`.
  - Modify `greville` function for the maximum multiplicity case to avoid
having the discontinuity points of the spline.

* Changes in modules `psydac.feec.derivatives` and
`psydac.feec.global_projectors` to take in account that the ghost
regions have size `shifts * pads` (previously was `pads`)

* Changes in `psydac.fem.splines`: repercussion of the changes in
`psydac.core.bsplines`, adding `multiplicity` as argument to several
functions. Correct the computation of the `multiplicity` from the knot
sequence.

* Add tests in `feec/tests/test_differentiation_matrices.py` and
`feec/tests/test_global_projectors.py` with higher multiplicity.

* Add test in `api/tests/tests_api_feec_1d.py,
api/tests/tests_api_feec_2d.py, api/tests/tests_api_feec_3d.py` with
higher multiplicity.

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.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
 
 
 def test_maxwell_1d_dirichlet_strong():

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,91 @@ 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_multiplicity():
+
+    namespace = run_maxwell_2d_TE(
+        use_spline_mapping = False,
+        eps      = 0.5,
+        ncells   = 10,
+        degree   = 5,
+        periodic = False,
+        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-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
+    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   = 30,
+        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
+    )
+
+    TOL = 1e-6
+    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
+    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 = 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
+    assert abs(namespace['error_l2_Bz'] - ref['error_l2_Bz']) / ref['error_l2_Bz'] <= TOL
+
 
 def test_maxwell_2d_dirichlet():
 

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)
@@ -356,6 +356,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.24749763720543216) < 1e-9
 
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
@@ -368,4 +408,8 @@ def teardown_module():
 def teardown_function():
     from sympy.core import cache
     cache.clear_cache()
+
+if __name__ == '__main__' :
+    test_maxwell_3d_2_mult()
+
 

psydac/core/bsplines.py

@@ -301,7 +301,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.
@@ -326,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
     -------
@@ -342,19 +346,20 @@ 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:
         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
 
 #==============================================================================
-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.
@@ -375,6 +380,10 @@ def histopolation_matrix(knots, degree, periodic, normalization, xgrid, check_bo
 
     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``.
@@ -418,23 +427,23 @@ 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"
 
     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')
-
-    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
 
 #==============================================================================
@@ -472,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.
 
@@ -490,6 +499,10 @@ def greville(knots, degree, periodic, 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
     -------
@@ -499,9 +512,10 @@ def greville(knots, degree, periodic, out=None):
     """
     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)
-    greville_p(knots, degree, periodic, out)
+    multiplicity = int(multiplicity)
+    greville_p(knots, degree, periodic, out, multiplicity)
     return out
 
 #===============================================================================
@@ -560,9 +574,12 @@ 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.
-    Otherwise, knot sequence is clamped (i.e. endpoints are repeated p times).
+    
+    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).
 
     Parameters
     ----------