diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index bdccc0e87..1ab23ec91 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,33 @@ 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] + 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], [1]) + V_cart = StencilVectorSpace(cart_decomp) + M = StencilMatrix(V_cart, V_cart) if cell == 'I': # interpolation case @@ -143,6 +159,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 + 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]) + + matrixcells += [M.copy()] + elif cell == 'H': # histopolation case if quad_x[j] is None: @@ -159,11 +192,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 + 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 +232,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 +304,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..07ed36d77 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -3,12 +3,12 @@ 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 mpi4py import MPI import numpy as np @@ -76,6 +76,29 @@ 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 + # TODO: test takes too long in 3D + #-------------------------- + # check BlockLinearOperator + #-------------------------- + # build the solver from the LinearOperator + # imat_kronecker_P0 = P0.imat_kronecker + # imat_kronecker_P1 = P1.imat_kronecker + # I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + # I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + + # # build the rhs + # P0.func(fun1) + # P1.func(D1fun1, D2fun1, D3fun1) + + # # solve and compare + # u0vec = u0.coeffs + # u0vec_imat = I0inv.solve(P0._rhs) + # assert np.allclose(u0vec.toarray(), u0vec_imat.toarray(), atol=1e-5) + + # u1vec = u1.coeffs + # u1vec_imat = I1inv.solve(P1._rhs) + # assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) + #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @@ -153,6 +176,29 @@ 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 + # TODO: test takes too long in 3D + #-------------------------- + # check BlockLinearOperator + #-------------------------- + # build the solver from the LinearOperator + # imat_kronecker_P1 = P1.imat_kronecker + # imat_kronecker_P2 = P2.imat_kronecker + # I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + # I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) + + # # build the rhs + # P1.func(fun1, fun2, fun3) + # P2.func(cf1, cf2, cf3) + + # # solve and compare + # u1vec = u1.coeffs + # u1vec_imat = I1inv.solve(P1._rhs) + # assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) + + # u2vec = u2.coeffs + # u2vec_imat = I2inv.solve(P2._rhs) + # assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5) + #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @@ -221,6 +267,30 @@ 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 + # TODO: test takes too long in 3D + #-------------------------- + # check BlockLinearOperator + #-------------------------- + # build the solver from the LinearOperator + # imat_kronecker_P2 = P2.imat_kronecker + # imat_kronecker_P3 = P3.imat_kronecker + # I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) + # I3inv = inverse(imat_kronecker_P3, 'gmres', verbose=True) + + # # build the rhs + # P2.func(fun1, fun2, fun3) + # P3.func(difun) + + # # solve and compare + # u2vec = u2.coeffs + # u2vec_imat = I2inv.solve(P2._rhs) + # assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5) + + # u3vec = u3.coeffs + # u3vec_imat = I3inv.solve(P3._rhs) + # assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-5) + +@pytest.mark.parallel @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @pytest.mark.parametrize('p', [2,3]) @@ -277,6 +347,29 @@ 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 + #-------------------------- + # build the solver from the LinearOperator + imat_kronecker_P0 = P0.imat_kronecker + imat_kronecker_P1 = P1.imat_kronecker + I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + + # build the rhs + P0.func(fun1) + P1.func(D1fun1, D2fun1) + + # solve and compare + u0vec = u0.coeffs + u0vec_imat = I0inv.solve(P0._rhs) + assert np.allclose(u0vec.toarray(), u0vec_imat.toarray(), atol=1e-5) + + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) + +@pytest.mark.parallel @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @pytest.mark.parametrize('p', [2,3]) @@ -333,6 +426,29 @@ 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 + #-------------------------- + # build the solver from the LinearOperator + imat_kronecker_P0 = P0.imat_kronecker + imat_kronecker_P1 = P1.imat_kronecker + I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + + # build the rhs + P0.func(fun1) + P1.func(D2fun1, D1fun1) + + # solve and compare + u0vec = u0.coeffs + u0vec_imat = I0inv.solve(P0._rhs) + assert np.allclose(u0vec.toarray(), u0vec_imat.toarray(), atol=1e-5) + + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) + +@pytest.mark.parallel @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @pytest.mark.parametrize('p', [2,3]) @@ -396,6 +512,29 @@ 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 + #-------------------------- + # build the solver from the LinearOperator + imat_kronecker_P2 = P2.imat_kronecker + imat_kronecker_P3 = P3.imat_kronecker + I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) + I3inv = inverse(imat_kronecker_P3, 'gmres', verbose=True) + + # build the rhs + P2.func(fun1, fun2) + P3.func(difun) + + # solve and compare + u2vec = u2.coeffs + u2vec_imat = I2inv.solve(P2._rhs) + assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5) + + u3vec = u3.coeffs + u3vec_imat = I3inv.solve(P3._rhs) + assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-5) + +@pytest.mark.parallel @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @pytest.mark.parametrize('p', [2,3]) @@ -451,14 +590,36 @@ 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 + #-------------------------- + # build the solver from the LinearOperator + imat_kronecker_P1 = P1.imat_kronecker + imat_kronecker_P2 = P2.imat_kronecker + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) + + # build the rhs + P1.func(fun1, fun2) + P2.func(difun) + + # solve and compare + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) + + u2vec = u2.coeffs + u2vec_imat = I2inv.solve(P2._rhs) + assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5) + @pytest.mark.parametrize('Nel', [16, 20]) @pytest.mark.parametrize('Nq', [5]) @pytest.mark.parametrize('p', [2,3]) @@ -509,6 +670,29 @@ 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 + #-------------------------- + # build the solver from the LinearOperator + imat_kronecker_P0 = P0.imat_kronecker + imat_kronecker_P1 = P1.imat_kronecker + I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + + # build the rhs + P0.func(fun1) + P1.func(Dfun1) + + # solve and compare + u0vec = u0.coeffs + u0vec_imat = I0inv.solve(P0._rhs) + assert np.allclose(u0vec.toarray(), u0vec_imat.toarray(), atol=1e-5) + + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) + #============================================================================== if __name__ == '__main__': 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