Expose inter-/histopolation matrices in GlobalProjector, add unit tests

psydac/feec/global_projectors.py

@@ -2,16 +2,17 @@
 
 import numpy as np
 
-from psydac.linalg.kron           import KroneckerLinearSolver
-from psydac.linalg.stencil        import StencilVector
-from psydac.linalg.block          import BlockLinearOperator, BlockVector
+from psydac.linalg.kron           import KroneckerLinearSolver, KroneckerStencilMatrix
+from psydac.linalg.stencil        import StencilMatrix, StencilVectorSpace
+from psydac.linalg.block          import BlockLinearOperator
 from psydac.core.bsplines         import quadrature_grid
 from psydac.utilities.quadratures import gauss_legendre
 from psydac.fem.basic             import FemField
 
 from psydac.fem.tensor import TensorFemSpace
 from psydac.fem.vector import VectorFemSpace
 
+from psydac.ddm.cart import DomainDecomposition, CartDecomposition
 from psydac.utilities.utils import roll_edges
 
 from abc import ABCMeta, abstractmethod
@@ -120,18 +121,34 @@ def __init__(self, space, nquads = None):
         self._grid_x = []
         self._grid_w = []
         solverblocks = []
+        matrixblocks = []
+        blocks = [] # for BlockLinearOperator
         for i,block in enumerate(structure):
             # do for each block (i.e. each TensorFemSpace):
-             
+
             block_x = []
             block_w = []
             solvercells = []
+            matrixcells = []
+            blocks += [[]] 
             for j, cell in enumerate(block):
                 # for each direction in the tensor space (i.e. each SplineSpace):
 
                 V = tensorspaces[i].spaces[j]
                 s = tensorspaces[i].vector_space.starts[j]
                 e = tensorspaces[i].vector_space.ends[j]
+                p = tensorspaces[i].vector_space.pads[j]
+                n = tensorspaces[i].vector_space.npts[j]
+                m = tensorspaces[i].multiplicity[j]
+                periodic = tensorspaces[i].vector_space.periods[j]
+                ncells = tensorspaces[i].ncells[j]
+                blocks[-1] += [None] # fill blocks with None, fill the diagonals later
+
+                # create a distributed matrix for the current 1d SplineSpace
+                domain_decomp = DomainDecomposition([ncells], [periodic])
+                cart_decomp = CartDecomposition(domain_decomp, [n], [[s]], [[e]], [p], [m])
+                V_cart = StencilVectorSpace(cart_decomp)
+                M = StencilMatrix(V_cart, V_cart)
 
                 if cell == 'I':
                     # interpolation case
@@ -143,6 +160,23 @@ def __init__(self, space, nquads = None):
                     local_x = local_intp_x[:, np.newaxis]
                     local_w = np.ones_like(local_x)
                     solvercells += [V._interpolator]
+
+                    # make 1D collocation matrix in stencil format
+                    row_indices, col_indices = np.nonzero(V.imat)
+
+                    for row_i, col_i in zip(row_indices, col_indices):
+
+                        # only consider row indices on process
+                        if row_i in range(V_cart.starts[0], V_cart.ends[0] + 1):
+                            row_i_loc = row_i - s
+
+                            M._data[row_i_loc + m*p, (col_i + p - row_i)%V.imat.shape[1]] = V.imat[row_i, col_i]
+
+                    # check if stencil matrix was built correctly
+                    assert np.allclose(M.toarray()[s:e + 1], V.imat[s:e + 1])
+                    # TODO Fix toarray() for multiplicity m > 1
+                    matrixcells += [M.copy()]
+
                 elif cell == 'H':
                     # histopolation case
                     if quad_x[j] is None:
@@ -159,11 +193,37 @@ def __init__(self, space, nquads = None):
 
                     local_x, local_w = quad_x[j], quad_w[j]
                     solvercells += [V._histopolator]
+
+                    # make 1D collocation matrix in stencil format
+                    row_indices, col_indices = np.nonzero(V.hmat)
+
+                    for row_i, col_i in zip(row_indices, col_indices):
+
+                        # only consider row indices on process
+                        if row_i in range(V_cart.starts[0], V_cart.ends[0] + 1):
+                            row_i_loc = row_i - s
+
+                            M._data[row_i_loc + m*p, (col_i + p - row_i)%V.hmat.shape[1]] = V.hmat[row_i, col_i]
+
+                    # check if stencil matrix was built correctly
+                    assert np.allclose(M.toarray()[s:e + 1], V.hmat[s:e + 1])
+
+                    matrixcells += [M.copy()]
+
                 else:
                     raise NotImplementedError('Invalid entry in structure array.')
 
                 block_x += [local_x]
                 block_w += [local_w]
+
+            # build Kronecker out of single directions    
+            if isinstance(self.space, TensorFemSpace):
+                matrixblocks += [KroneckerStencilMatrix(self.space.vector_space, self.space.vector_space, *matrixcells)]
+            else:
+                matrixblocks += [KroneckerStencilMatrix(self.space.vector_space[i], self.space.vector_space[i], *matrixcells)]
+
+            # fill the diagonals for BlockLinearOperator
+            blocks[i][i] = matrixblocks[-1]
 
             # finish block, build solvers, get dataslice to project to
             self._grid_x += [block_x]
@@ -173,6 +233,13 @@ def __init__(self, space, nquads = None):
 
             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]
+
+        # build final Inter-/Histopolation matrix (distributed)        
+        if isinstance(self.space, TensorFemSpace):
+            self._imat_kronecker = matrixblocks[0]
+        else:
+            self._imat_kronecker = BlockLinearOperator(self.space.vector_space, self.space.vector_space, 
+                                                       blocks=blocks)
 
         # finish arguments and create a lambda
         args = (*intp_x, *quad_x, *quad_w, *dofs)
@@ -238,6 +305,13 @@ def solver(self):
         """
         return self._solver
 
+    @property
+    def imat_kronecker(self):
+        """
+        Inter-/Histopolation matrix in distributed format.
+        """
+        return self._imat_kronecker
+
     @abstractmethod
     def _structure(self, dim):
         """

psydac/feec/tests/test_commuting_projections.py

@@ -3,19 +3,20 @@
 from psydac.feec.global_projectors import Projector_H1, Projector_L2, Projector_Hcurl, Projector_Hdiv
 from psydac.fem.tensor       import TensorFemSpace, SplineSpace
 from psydac.fem.vector       import VectorFemSpace
-from psydac.linalg.block     import BlockVector
 from psydac.core.bsplines    import make_knots
 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.ddm.cart         import DomainDecomposition
+from psydac.linalg.solvers   import inverse
+from psydac.linalg.basic     import IdentityOperator
 
 from mpi4py import MPI
 import numpy as np
 import pytest
 
 #==============================================================================
-# Run test
+# 3D tests
 #==============================================================================
 @pytest.mark.parametrize('Nel', [8, 12])
 @pytest.mark.parametrize('Nq', [5])
@@ -76,7 +77,21 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc, m):
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
 
-#==============================================================================
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+    Id_0 = IdentityOperator(H1.vector_space)
+    Err_0 = P0.solver @ P0.imat_kronecker - Id_0
+    e0 = Err_0 @ u0.coeffs  # random vector could be used as well
+    norm2_e0 = np.sqrt(e0.dot(e0))
+    assert norm2_e0 < 1e-12
+
+    Id_1 = IdentityOperator(Hcurl.vector_space)
+    Err_1 = P1.solver @ P1.imat_kronecker - Id_1
+    e1 = Err_1 @ u1.coeffs  # random vector could be used as well
+    norm2_e1 = np.sqrt(e1.dot(e1))
+    assert norm2_e1 < 1e-12
+
 @pytest.mark.parametrize('Nel', [8, 12])
 @pytest.mark.parametrize('Nq', [8])
 @pytest.mark.parametrize('p', [2,3])
@@ -153,7 +168,22 @@ def test_3d_commuting_pro_2(Nel, Nq, p, bc, m):
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
 
-#==============================================================================
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+
+    Id_1 = IdentityOperator(Hcurl.vector_space)
+    Err_1 = P1.solver @ P1.imat_kronecker - Id_1
+    e1 = Err_1 @ u1.coeffs  # random vector could be used as well
+    norm2_e1 = np.sqrt(e1.dot(e1))
+    assert norm2_e1 < 1e-12
+
+    Id_2 = IdentityOperator(Hdiv.vector_space)
+    Err_2 = P2.solver @ P2.imat_kronecker - Id_2
+    e2 = Err_2 @ u2.coeffs  # random vector could be used as well
+    norm2_e2 = np.sqrt(e2.dot(e2))
+    assert norm2_e2 < 1e-12
+
 @pytest.mark.parametrize('Nel', [8, 12])
 @pytest.mark.parametrize('Nq', [8])
 @pytest.mark.parametrize('p', [2,3])
@@ -221,6 +251,26 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc, m):
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
 
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+
+    Id_2 = IdentityOperator(Hdiv.vector_space)
+    Err_2 = P2.solver @ P2.imat_kronecker - Id_2
+    e2 = Err_2 @ u2.coeffs  # random vector could be used as well
+    norm2_e2 = np.sqrt(e2.dot(e2))
+    assert norm2_e2 < 1e-12
+
+    Id_3 = IdentityOperator(L2.vector_space)
+    Err_3 = P3.solver @ P3.imat_kronecker - Id_3
+    e3 = Err_3 @ u3.coeffs  # random vector could be used as well
+    norm2_e3 = np.sqrt(e3.dot(e3))
+    assert norm2_e3 < 1e-12
+
+#==============================================================================
+# 2D tests
+#==============================================================================
+@pytest.mark.parallel
 @pytest.mark.parametrize('Nel', [8, 12])
 @pytest.mark.parametrize('Nq', [5])
 @pytest.mark.parametrize('p', [2,3])
@@ -277,6 +327,23 @@ def test_2d_commuting_pro_1(Nel, Nq, p, bc, m):
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
 
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+
+    Id_0 = IdentityOperator(H1.vector_space)
+    Err_0 = P0.solver @ P0.imat_kronecker - Id_0
+    e0 = Err_0 @ u0.coeffs  # random vector could be used as well
+    norm2_e0 = np.sqrt(e0.dot(e0))
+    assert norm2_e0 < 1e-12
+
+    Id_1 = IdentityOperator(Hcurl.vector_space)
+    Err_1 = P1.solver @ P1.imat_kronecker - Id_1
+    e1 = Err_1 @ u1.coeffs  # random vector could be used as well
+    norm2_e1 = np.sqrt(e1.dot(e1))
+    assert norm2_e1 < 1e-12
+
+@pytest.mark.parallel
 @pytest.mark.parametrize('Nel', [8, 12])
 @pytest.mark.parametrize('Nq', [5])
 @pytest.mark.parametrize('p', [2,3])
@@ -333,6 +400,23 @@ def test_2d_commuting_pro_2(Nel, Nq, p, bc, m):
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
 
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+
+    Id_0 = IdentityOperator(H1.vector_space)
+    Err_0 = P0.solver @ P0.imat_kronecker - Id_0
+    e0 = Err_0 @ u0.coeffs  # random vector could be used as well
+    norm2_e0 = np.sqrt(e0.dot(e0))
+    assert norm2_e0 < 1e-12
+
+    Id_1 = IdentityOperator(Hdiv.vector_space)
+    Err_1 = P1.solver @ P1.imat_kronecker - Id_1
+    e1 = Err_1 @ u1.coeffs  # random vector could be used as well
+    norm2_e1 = np.sqrt(e1.dot(e1))
+    assert norm2_e0 < 1e-12
+
+@pytest.mark.parallel
 @pytest.mark.parametrize('Nel', [8, 12])
 @pytest.mark.parametrize('Nq', [8])
 @pytest.mark.parametrize('p', [2,3])
@@ -396,6 +480,23 @@ def test_2d_commuting_pro_3(Nel, Nq, p, bc, m):
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
 
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+
+    Id_2 = IdentityOperator(Hdiv.vector_space)
+    Err_2 = P2.solver @ P2.imat_kronecker - Id_2
+    e2 = Err_2 @ u2.coeffs
+    norm2_e2 = np.sqrt(e2.dot(e2))
+    assert norm2_e2 < 1e-12
+
+    Id_3 = IdentityOperator(L2.vector_space)
+    Err_3 = P3.solver @ P3.imat_kronecker - Id_3
+    e3 = Err_3 @ u3.coeffs
+    norm2_e3 = np.sqrt(e3.dot(e3))
+    assert norm2_e3 < 1e-12
+
+@pytest.mark.parallel
 @pytest.mark.parametrize('Nel', [8, 12])
 @pytest.mark.parametrize('Nq', [8])
 @pytest.mark.parametrize('p', [2,3])
@@ -451,14 +552,33 @@ def test_2d_commuting_pro_4(Nel, Nq, p, bc, m):
     # Projections and discrete derivatives
     #-------------------------------------
 
-    u2        = P1((fun1, fun2))
-    u3        = P2(difun)
-    Dfun_h    = curl(u2)
-    Dfun_proj = u3
+    u1        = P1((fun1, fun2))
+    u2        = P2(difun)
+    Dfun_h    = curl(u1)
+    Dfun_proj = u2
 
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
 
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+
+    Id_1 = IdentityOperator(Hcurl.vector_space)
+    Err_1 = P1.solver @ P1.imat_kronecker - Id_1
+    e1 = Err_1 @ u1.coeffs
+    norm2_e1 = np.sqrt(e1.dot(e1))
+    assert norm2_e1 < 1e-12
+
+    Id_2 = IdentityOperator(L2.vector_space)
+    Err_2 = P2.solver @ P2.imat_kronecker - Id_2
+    e2 = Err_2 @ u2.coeffs
+    norm2_e2 = np.sqrt(e2.dot(e2))
+    assert norm2_e2 < 1e-12
+
+#==============================================================================
+# 1D tests
+#==============================================================================
 @pytest.mark.parametrize('Nel', [16, 20])
 @pytest.mark.parametrize('Nq', [5])
 @pytest.mark.parametrize('p', [2,3])
@@ -509,6 +629,23 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc, m):
 
     error = abs((Dfun_proj.coeffs-Dfun_h.coeffs).toarray()).max()
     assert error < 1e-9
+
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+
+    Id_0 = IdentityOperator(H1.vector_space)
+    Err_0 = P0.solver @ P0.imat_kronecker - Id_0
+    e0 = Err_0 @ u0.coeffs
+    norm2_e0 = np.sqrt(e0.dot(e0))
+    assert norm2_e0 < 1e-12
+
+    Id_1 = IdentityOperator(L2.vector_space)
+    Err_1 = P1.solver @ P1.imat_kronecker - Id_1
+    e1 = Err_1 @ u1.coeffs
+    norm2_e1 = np.sqrt(e1.dot(e1))
+    assert norm2_e1 < 1e-12
+
 #==============================================================================
 if __name__ == '__main__':
 
@@ -526,4 +663,3 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc, m):
     test_2d_commuting_pro_3(Nel, Nq, p, bc, m)
     test_2d_commuting_pro_4(Nel, Nq, p, bc, m)
     test_1d_commuting_pro_1(Nel, Nq, p, bc, m)
-

psydac/fem/splines.py

@@ -193,7 +193,7 @@ def init_interpolation( self, dtype=float ):
             for i,j in zip( *cmat.nonzero() ):
                 bmat[u+l+i-j,j] = cmat[i,j]
             self._interpolator = BandedSolver( u, l, bmat )
-            self.imat = imat
+        self.imat = imat
 
         # Store flag
         self._interpolation_ready = True