diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index bdccc0e87..509995406 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -2,9 +2,9 @@ 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 @@ -12,6 +12,7 @@ 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): """ diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 025cfadd7..006439b01 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/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) - diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 11ebbc1f2..bb475ecb2 100644 --- a/psydac/fem/splines.py +++ b/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 diff --git a/psydac/linalg/kron.py b/psydac/linalg/kron.py index d68293c0c..e4abfdd33 100644 --- a/psydac/linalg/kron.py +++ b/psydac/linalg/kron.py @@ -94,19 +94,19 @@ def dot(self, x, out=None): starts = self._codomain.starts ends = self._codomain.ends pads = self._codomain.pads + shifts = self._codomain.shifts mats = self.mats - nrows = tuple(e-s+1 for s,e in zip(starts, ends)) pnrows = tuple(2*p+1 for p in pads) - + for ii in np.ndindex(*nrows): v = 0. - xx = tuple(i+p for i,p in zip(ii, pads)) + xx = tuple(i+p*s for i,p,s in zip(ii, pads, shifts)) for jj in np.ndindex(*pnrows): i_mats = [mat._data[s, j] for s,j,mat in zip(xx, jj, mats)] - ii_jj = tuple(i+j for i,j in zip(ii, jj)) + ii_jj = tuple(i+j+(s-1)*p for i,j,p,s in zip(ii, jj, pads, shifts)) v += x._data[ii_jj] * np.prod(i_mats) out._data[xx] = v