Merge branch 'devel' into expose-inter-histopolation-matrices

docs/source/modules/linalg/block.rst

@@ -4,7 +4,6 @@ linalg.block
    * :ref:`BlockVectorSpace <blockvectorspace>`
    * :ref:`BlockVector <blockvector>`
    * :ref:`BlockLinearOperator <blocklinearoperator>`
-   * :ref:`BlockDiagonalSolver <blockdiagonalsolver>`  
 
 .. inheritance-diagram:: psydac.linalg.block
 
@@ -31,11 +30,3 @@ BlockLinearOperator
 
 .. autoclass:: psydac.linalg.block.BlockLinearOperator
    :members:
-
-.. _blockdiagonalsolver:
-
-BlockDiagonalSolver
--------------------
-
-.. autoclass:: psydac.linalg.block.BlockDiagonalSolver
-   :members:

docs/source/modules/linalg/direct_solvers.rst

@@ -1,20 +1,11 @@
 linalg.direct_solvers
 =====================
 
-   * :ref:`DirectSolver <directsolver>`
    * :ref:`BandedSolver <bandedsolver>`
    * :ref:`SparseSolver <sparsesolver>`
 
 .. inheritance-diagram:: psydac.linalg.direct_solvers
 
-.. _directsolver:
-
-DirectSolver
-------------
-
-.. autoclass:: psydac.linalg.direct_solvers.DirectSolver
-   :members:
-
 .. _bandedsolver:
 
 BandedSolver

psydac/api/ast/evaluation.py

@@ -173,8 +173,6 @@ def _initialize(self):
         header = None
         if self.backend['name'] == 'pyccel':
             decorators = {'types': build_pyccel_types_decorator(func_args)}
-        elif self.backend['name'] == 'numba':
-            decorators = {'jit':[]}
         elif self.backend['name'] == 'pythran':
             header = build_pythran_types_header(self.name, func_args)
 
@@ -410,8 +408,6 @@ def _initialize(self):
         header = None
         if self.backend['name'] == 'pyccel':
             decorators = {'types': build_pyccel_types_decorator(func_args)}
-        elif self.backend['name'] == 'numba':
-            decorators = {'jit':[]}
         elif self.backend['name'] == 'pythran':
             header = build_pythran_types_header(self.name, func_args)
 

psydac/api/ast/expr.py

@@ -434,8 +434,6 @@ def _initialize(self):
         header = None
         if self.backend['name'] == 'pyccel':
             decorators = {'types': build_pyccel_types_decorator(func_args)}
-        elif self.backend['name'] == 'numba':
-            decorators = {'jit':[]}
         elif self.backend['name'] == 'pythran':
             header = build_pythran_types_header(self.name, func_args)
 

psydac/api/ast/glt.py

@@ -562,8 +562,6 @@ def _initialize(self, **kwargs):
         header = None
         if self.backend['name'] == 'pyccel':
             decorators = {'types': build_pyccel_types_decorator(func_args)}
-        elif self.backend['name'] == 'numba':
-            decorators = {'jit':[]}
         elif self.backend['name'] == 'pythran':
             header = build_pythran_types_header(self.name, func_args)
 

psydac/api/ast/linalg.py

@@ -345,8 +345,6 @@ def _initialize(self, ndim, block_shape, **kwargs):
             if backend['name'] == 'pyccel':
                 a = [String(str(i)) for i in build_pyccel_types_decorator(func_args)]
                 decorators = {'types': Function('types')(*a)}
-            elif backend['name'] == 'numba':
-                decorators = {'njit': Function('njit')(ValuedArgument(Symbol('fastmath'), backend['fastmath']))}
             elif backend['name'] == 'pythran':
                 header = build_pythran_types_header(name, func_args)
 
@@ -390,9 +388,6 @@ def _generate_code(self, backend=None):
             if backend and backend['name'] == 'pyccel':
                 imports  = 'from pyccel.decorators import types\n'
                 imports += 'from numpy import shape'
-            elif backend and backend['name'] == 'numba':
-                imports  = 'from numba import njit\n'
-                imports += 'from numpy import shape'
             else:
                 imports = 'from numpy import shape'
 
@@ -494,8 +489,6 @@ def _initialize(self):
 
         if self.backend['name'] == 'pyccel':
             decorators = {'types': build_pyccel_types_decorator(func_args), 'external':[]}
-        elif self.backend['name'] == 'numba':
-            decorators = {'jit':[]}
         elif self.backend['name'] == 'pythran':
             header = build_pythran_types_header(self.name, func_args)
 

psydac/api/ast/parser.py

@@ -562,19 +562,13 @@ def _visit_DefNode(self, expr, **kwargs):
         if self.backend['name'] == 'pyccel':
             a = [String(str(i)) for i in build_pyccel_types_decorator(arguments)]
             decorators = {'types': Function('types')(*a)}
-        elif self.backend['name'] == 'numba':
-            decorators = {'njit': Function('njit')(ValuedArgument(Symbol('fastmath'), self.backend['fastmath']))}
         elif self.backend['name'] == 'pythran':
             header = build_pythran_types_header(name, arguments)
         else:
             decorators = {}
 
-        if self.backend['name'] == 'numba':
-            func = FunctionDef(name, arguments, results, body, decorators=decorators)
-            stmts = CodeBlock([*imports , func])
-        else:
-            func = FunctionDef(name, arguments, results, body, imports=imports, decorators=decorators)
-            stmts = func
+        func = FunctionDef(name, arguments, results, body, imports=imports, decorators=decorators)
+        stmts = func
 
         self.functions[name] = func
         return stmts

psydac/api/ast/utilities.py

@@ -331,14 +331,14 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad,
         orders = [*get_index(atom).values()]
         args   = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
         inits += [Assign(test_fun, Mul(*args))]
-        # ...
+        # ...
 
-        # ...
+        # ...
         args     = [IndexedBase(field_name)[idxs], test_fun]
         val_name = SymbolicExpr(atom).name + '_values'
         val      = IndexedBase(val_name)[indices_quad]
         updates += [AugAssign(val,'+',Mul(*args))]
-        # ...
+        # ...
 
     return inits, updates, map_stmts, new_atoms
 
@@ -386,23 +386,23 @@ def compute_atoms_expr_mapping(atomic_exprs, indices_quad,
         element = get_atom_logical_derivatives(atom)
         element_name = 'coeff_' + SymbolicExpr(element).name
 
-        # ...
+        # ...
         test_fun = atom.subs(element, test_function)
         test_fun = SymbolicExpr(test_fun)
-        # ...
+        # ...
 
-        # ...
+        # ...
         orders = [*get_index_logical_derivatives(atom).values()]
         args   = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
         inits += [Assign(test_fun, Mul(*args))]
-        # ...
+        # ...
 
-        # ...
+        # ...
         val_name = SymbolicExpr(atom).name + '_values'
         val      = IndexedBase(val_name)[indices_quad]
         expr     = IndexedBase(element_name)[idxs] * test_fun
         updates += [AugAssign(val, '+', expr)]
-        # ...
+        # ...
 
     return inits, updates
 

psydac/api/basic.py

@@ -229,11 +229,10 @@ def _initialize_folder(self, folder=None):
     def _generate_code(self):
         # ... generate code that can be pyccelized
         imports = ''
+
         if self.backend['name'] == 'pyccel':
             imports = "from pyccel.decorators import types"
             imports += ", template \n@template(name='T', types=['float', 'complex']) "
-        elif self.backend['name'] == 'numba':
-            imports = 'from numba import njit'
 
         ast = self.ast
         expr = parse(ast.expr, settings={'dim': ast.dim, 'nderiv': ast.nderiv, 'mapping':ast.mapping, 'target':ast.domain}, backend=self.backend)

psydac/api/discretization.py

@@ -81,18 +81,42 @@ def change_dtype(V, dtype):
 
     return V
 
-#==============================================================================
-def discretize_derham(derham, domain_h, *args, **kwargs):
+#==============================================================================           
+def discretize_derham(derham, domain_h, get_H1vec_space = False, *args, **kwargs):
+    """
+    Create a discrete De Rham sequence by creating the spaces and then initiating DiscreteDerham object.
+
+    Parameters
+    ----------
+
+    derham : sympde.topology.space.Derham
+        The symbolic Derham sequence
+
+    domain_h   : Geometry
+        Discrete domain where the spaces will be discretized
+        
+    get_H1vec_space : Bool
+        True to also get the "Hvec" space discretizing (H1)^n vector fields
+        
+    **kwargs : list
+        optional parameters for the space discretization
+    """
 
     ldim    = derham.shape
     mapping = domain_h.domain.mapping # NOTE: assuming single-patch domain!
-
     bases  = ['B'] + ldim * ['M']
     spaces = [discretize_space(V, domain_h, basis=basis, **kwargs) \
             for V, basis in zip(derham.spaces, bases)]
 
-    return DiscreteDerham(mapping, *spaces)
+    if get_H1vec_space:
+        X = VectorFunctionSpace('X', domain_h.domain, kind='h1')
+        V0h = spaces[0]
+        Xh  = VectorFemSpace(*([V0h]*ldim))
+        Xh.symbolic_space = X
+        #We still need to specify the symbolic space because of "_recursive_element_of" not implemented in sympde
+        spaces.append(Xh)
 
+    return DiscreteDerham(mapping, *spaces)
 #==============================================================================
 def reduce_space_degrees(V, Vh, *, basis='B', sequence='DR'):
     """

psydac/api/feec.py

@@ -4,25 +4,46 @@
 from psydac.feec.derivatives       import Derivative_1D, Gradient_2D, Gradient_3D
 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, Projector_Hcurl
+from psydac.feec.global_projectors import Projector_H1, Projector_Hcurl, Projector_H1vec
 from psydac.feec.global_projectors import Projector_Hdiv, Projector_L2
 from psydac.feec.pull_push         import pull_1d_h1, pull_1d_l2
-from psydac.feec.pull_push         import pull_2d_h1, pull_2d_hcurl, pull_2d_hdiv, pull_2d_l2
-from psydac.feec.pull_push         import pull_3d_h1, pull_3d_hcurl, pull_3d_hdiv, pull_3d_l2
+from psydac.feec.pull_push         import pull_2d_h1, pull_2d_hcurl, pull_2d_hdiv, pull_2d_l2, pull_2d_h1vec
+from psydac.feec.pull_push         import pull_3d_h1, pull_3d_hcurl, pull_3d_hdiv, pull_3d_l2, pull_3d_h1vec
+from psydac.fem.vector             import VectorFemSpace
+
 
 __all__ = ('DiscreteDerham',)
 
 #==============================================================================
 class DiscreteDerham(BasicDiscrete):
     """ Represent the discrete De Rham sequence.
+    Should be initialized via discretize_derham function in api.discretization.py
+    
+    Parameters
+    ----------
+
+    mapping : Mapping
+        The mapping from the logical space to the physical space of the discrete De Rham.
+        
+    *spaces : list of FemSpace
+        The discrete spaces of the De Rham sequence
     """
     def __init__(self, mapping, *spaces):
 
         assert (mapping is None) or isinstance(mapping, Mapping)
+
+        self.has_vec = isinstance(spaces[-1], VectorFemSpace)
+
+        if self.has_vec : 
+            dim          = len(spaces) - 2
+            self._spaces = spaces[:-1]
+            self._H1vec    = spaces[-1]
+
+        else :
+            dim           = len(spaces) - 1
+            self._spaces  = spaces
 
-        dim           = len(spaces) - 1
         self._dim     = dim
-        self._spaces  = spaces
         self._mapping = mapping
         self._callable_mapping = mapping.get_callable_mapping() if mapping else None
 
@@ -83,6 +104,11 @@ def V2(self):
     def V3(self):
         return self._spaces[3]
 
+    @property
+    def H1vec(self):
+        assert self.has_vec
+        return self._H1vec
+
     @property
     def spaces(self):
         return self._spaces
@@ -130,25 +156,47 @@ def projectors(self, *, kind='global', nquads=None):
             else:
                 raise TypeError('projector of space type {} is not available'.format(kind))
 
+            if self.has_vec : 
+                Pvec = Projector_H1vec(self.H1vec, nquads)
+
             if self.mapping:
                 P0_m = lambda f: P0(pull_2d_h1(f, self.callable_mapping))
                 P2_m = lambda f: P2(pull_2d_l2(f, self.callable_mapping))
                 if kind == 'hcurl':
                     P1_m = lambda f: P1(pull_2d_hcurl(f, self.callable_mapping))
                 elif kind == 'hdiv':
                     P1_m = lambda f: P1(pull_2d_hdiv(f, self.callable_mapping))
-                return P0_m, P1_m, P2_m
-            return P0, P1, P2
+                if self.has_vec : 
+                    Pvec_m = lambda f: Pvec(pull_2d_h1vec(f, self.callable_mapping))
+                    return P0_m, P1_m, P2_m, Pvec_m
+                else : 
+                    return P0_m, P1_m, P2_m
+
+            if self.has_vec :
+                return P0, P1, P2, Pvec
+            else : 
+                return P0, P1, P2
 
         elif self.dim == 3:
             P0 = Projector_H1   (self.V0)
             P1 = Projector_Hcurl(self.V1, nquads)
             P2 = Projector_Hdiv (self.V2, nquads)
             P3 = Projector_L2   (self.V3, nquads)
+            if self.has_vec : 
+                Pvec = Projector_H1vec(self.H1vec)
             if self.mapping:
                 P0_m = lambda f: P0(pull_3d_h1   (f, self.callable_mapping))
                 P1_m = lambda f: P1(pull_3d_hcurl(f, self.callable_mapping))
                 P2_m = lambda f: P2(pull_3d_hdiv (f, self.callable_mapping))
                 P3_m = lambda f: P3(pull_3d_l2   (f, self.callable_mapping))
-                return P0_m, P1_m, P2_m, P3_m
-            return P0, P1, P2, P3
+                if self.has_vec : 
+                    Pvec_m = lambda f: Pvec(pull_3d_h1vec(f, self.callable_mapping))
+                    return P0_m, P1_m, P2_m, P3_m, Pvec_m
+                else : 
+                    return P0_m, P1_m, P2_m, P3_m
+
+            if self.has_vec :
+                return P0, P1, P2, P3, Pvec
+            else : 
+                return P0, P1, P2, P3
+

psydac/api/fem.py

@@ -1539,7 +1539,7 @@ def assemble(self, **kwargs):
             v = kwargs[key]
             if isinstance(v, FemField):
                 if not v.coeffs.ghost_regions_in_sync:
-                        v.coeffs.update_ghost_regions()
+                    v.coeffs.update_ghost_regions()
                 if v.space.is_product:
                     coeffs = v.coeffs
                     if self._symbolic_space.is_broken:

psydac/api/glt.py

@@ -266,9 +266,6 @@ def _generate_code(self):
             code += '\nfrom pyccel.decorators import types'
             code += '\nfrom pyccel.decorators import external, external_call'
 
-        elif self.backend['name'] == 'numba':
-            code = 'from numba import jit'
-
         imports = '\n'.join(pycode(imp) for dep in self.dependencies for imp in dep.imports )
 
         code = '{code}\n{imports}'.format(code=code, imports=imports)

psydac/api/settings.py

@@ -32,10 +32,6 @@
                       'folder': '__pgpyccel__',
                        'tag':'pgpyccel',
                        'openmp':False}
-
-PSYDAC_BACKEND_NUMBA = {'name': 'numba','tag':'numba', 'fastmath':True, 'openmp':False}
-
-#PSYDAC_BACKEND_PYTHRAN = {'name':'pythran','tag':'pythran', 'openmp':False}
 # ...
 
 # List of all available backends for accelerating Python code
@@ -44,6 +40,4 @@
     'pyccel-gcc'  : PSYDAC_BACKEND_GPYCCEL,
     'pyccel-intel': PSYDAC_BACKEND_IPYCCEL,
     'pyccel-pgi'  : PSYDAC_BACKEND_PGPYCCEL,
-    'numba'       : PSYDAC_BACKEND_NUMBA,
- #   'pythran'     : PSYDAC_BACKEND_PYTHRAN,
 }