De Rham sequences w/ knot multiplicity > 1

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

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, points[0, 0], space.periodic)
                 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, xq, perio )
         glob_spans[iq] = span
 
         ders = basis_funs_all_ders( T, p, xq, span, nderiv )

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

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-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   = 10,
+        degree   = 5,
+        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 = 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
+    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():
 

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.9020674104318802) < 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()
+