From ba29e906917670822aa2aa5264b6ce58ea692933 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 12 Oct 2023 15:30:38 +0200 Subject: [PATCH 1/6] Expose inter-/histopolation matrices in GlobalProjector, add unit tests --- psydac/feec/global_projectors.py | 82 +++++++- .../feec/tests/test_commuting_projections.py | 197 +++++++++++++++++- psydac/fem/splines.py | 2 +- 3 files changed, 269 insertions(+), 12 deletions(-) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 773df5cf3..74b2abca3 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 BlockDiagonalSolver, BlockVector +from psydac.linalg.kron import KroneckerLinearSolver, KroneckerStencilMatrix +from psydac.linalg.stencil import StencilMatrix, StencilVectorSpace +from psydac.linalg.block import BlockDiagonalSolver, BlockLinearOperator from psydac.core.bsplines import quadrature_grid from psydac.utilities.quadratures import gauss_legendre from psydac.fem.basic import FemField @@ -12,6 +12,8 @@ from psydac.fem.tensor import TensorFemSpace from psydac.fem.vector import VectorFemSpace +from psydac.ddm.cart import DomainDecomposition, CartDecomposition + from abc import ABCMeta, abstractmethod __all__ = ('GlobalProjector', 'Projector_H1', 'Projector_Hcurl', 'Projector_Hdiv', 'Projector_L2', @@ -119,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 @@ -142,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: @@ -151,11 +185,37 @@ def __init__(self, space, nquads = None): quad_w[j] = global_quad_w[s:e+1] 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] @@ -165,6 +225,13 @@ def __init__(self, space, nquads = None): dataslice = tuple(slice(p, -p) for p in tensorspaces[i].vector_space.pads) 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) @@ -228,6 +295,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 f4cc6c2c3..2d612622e 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 GMRES as Inverse from mpi4py import MPI import numpy as np @@ -74,6 +74,28 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc): 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, verbose=True) + I1inv = Inverse(imat_kronecker_P1, 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-6) + + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-6) + #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @@ -149,6 +171,28 @@ def test_3d_commuting_pro_2(Nel, Nq, p, bc): 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, verbose=True) + I2inv = Inverse(imat_kronecker_P2, 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-6) + + u2vec = u2.coeffs + u2vec_imat = I2inv.solve(P2._rhs) + assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-6) + #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @@ -215,6 +259,28 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc): 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, verbose=True) + I3inv = Inverse(imat_kronecker_P3, 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-6) + + u3vec = u3.coeffs + u3vec_imat = I3inv.solve(P3._rhs) + assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-6) + @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @pytest.mark.parametrize('p', [2,3]) @@ -269,6 +335,28 @@ def test_2d_commuting_pro_1(Nel, Nq, p, bc): 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, verbose=True) + I1inv = Inverse(imat_kronecker_P1, 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()) + + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray()) + @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @pytest.mark.parametrize('p', [2,3]) @@ -323,6 +411,28 @@ def test_2d_commuting_pro_2(Nel, Nq, p, bc): 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, verbose=True) + I1inv = Inverse(imat_kronecker_P1, 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()) + + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray()) + @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @pytest.mark.parametrize('p', [2,3]) @@ -384,6 +494,28 @@ def test_2d_commuting_pro_3(Nel, Nq, p, bc): 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, verbose=True) + I3inv = Inverse(imat_kronecker_P3, 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-6) + + u3vec = u3.coeffs + u3vec_imat = I3inv.solve(P3._rhs) + assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-6) + @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @pytest.mark.parametrize('p', [2,3]) @@ -437,14 +569,36 @@ def test_2d_commuting_pro_4(Nel, Nq, p, bc): # 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, verbose=True) + I2inv = Inverse(imat_kronecker_P2, 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-6) + + u2vec = u2.coeffs + u2vec_imat = I2inv.solve(P2._rhs) + assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-6) + @pytest.mark.parametrize('Nel', [16, 20]) @pytest.mark.parametrize('Nq', [5]) @pytest.mark.parametrize('p', [2,3]) @@ -493,15 +647,44 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc): 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, verbose=True) + I1inv = Inverse(imat_kronecker_P1, 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()) + + u1vec = u1.coeffs + u1vec_imat = I1inv.solve(P1._rhs) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray()) + #============================================================================== if __name__ == '__main__': Nel = 8 Nq = 8 p = 3 - bc = True + bc = False + test_1d_commuting_pro_1(Nel, Nq, p, bc) + test_2d_commuting_pro_1(Nel, Nq, p, bc) + test_2d_commuting_pro_2(Nel, Nq, p, bc) + test_2d_commuting_pro_3(Nel, Nq, p, bc) + test_2d_commuting_pro_4(Nel, Nq, p, bc) test_3d_commuting_pro_1(Nel, Nq, p, bc) - test_3d_commuting_pro_2(Nel, Nq, p, bc) test_3d_commuting_pro_3(Nel, Nq, p, bc) + test_2d_commuting_pro_4(Nel, Nq, p, bc) + diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index e18d09c94..be08c44a5 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -194,7 +194,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 From 60592c988749ab84c90cbaa1ef10e94e646e43a3 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 12 Oct 2023 19:47:17 +0200 Subject: [PATCH 2/6] Set atol=1e-5 in tests --- .../feec/tests/test_commuting_projections.py | 39 ++++++++++--------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 2d612622e..4ad1f9591 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -90,11 +90,11 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc): # solve and compare u0vec = u0.coeffs u0vec_imat = I0inv.solve(P0._rhs) - assert np.allclose(u0vec.toarray(), u0vec_imat.toarray(), atol=1e-6) + 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-6) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) @@ -187,11 +187,11 @@ def test_3d_commuting_pro_2(Nel, Nq, p, bc): # solve and compare u1vec = u1.coeffs u1vec_imat = I1inv.solve(P1._rhs) - assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-6) + 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-6) + assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5) #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) @@ -275,11 +275,11 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc): # solve and compare u2vec = u2.coeffs u2vec_imat = I2inv.solve(P2._rhs) - assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-6) + 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-6) + assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-5) @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @@ -510,11 +510,11 @@ def test_2d_commuting_pro_3(Nel, Nq, p, bc): # solve and compare u2vec = u2.coeffs u2vec_imat = I2inv.solve(P2._rhs) - assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-6) + 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-6) + assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-5) @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @@ -593,11 +593,11 @@ def test_2d_commuting_pro_4(Nel, Nq, p, bc): # solve and compare u1vec = u1.coeffs u1vec_imat = I1inv.solve(P1._rhs) - assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-6) + 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-6) + assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5) @pytest.mark.parametrize('Nel', [16, 20]) @pytest.mark.parametrize('Nq', [5]) @@ -673,18 +673,19 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc): #============================================================================== if __name__ == '__main__': - Nel = 8 + Nel = 12 Nq = 8 p = 3 bc = False - test_1d_commuting_pro_1(Nel, Nq, p, bc) - test_2d_commuting_pro_1(Nel, Nq, p, bc) - test_2d_commuting_pro_2(Nel, Nq, p, bc) - test_2d_commuting_pro_3(Nel, Nq, p, bc) - test_2d_commuting_pro_4(Nel, Nq, p, bc) - test_3d_commuting_pro_1(Nel, Nq, p, bc) - test_3d_commuting_pro_3(Nel, Nq, p, bc) - test_2d_commuting_pro_4(Nel, Nq, p, bc) + # test_1d_commuting_pro_1(Nel, Nq, p, bc) + # test_2d_commuting_pro_1(Nel, Nq, p, bc) + # test_2d_commuting_pro_2(Nel, Nq, p, bc) + # test_2d_commuting_pro_3(Nel, Nq, p, bc) + # test_2d_commuting_pro_4(Nel, Nq, p, bc) + # test_3d_commuting_pro_1(Nel, Nq, p, bc) + test_3d_commuting_pro_2(Nel, Nq, p, bc) + # test_3d_commuting_pro_3(Nel, Nq, p, bc) + # test_2d_commuting_pro_4(Nel, Nq, p, bc) From b59195a1c6457fb88805fcace8e9236e159169b5 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 13 Oct 2023 08:06:33 +0200 Subject: [PATCH 3/6] Add forgotten atol, set pyccel backend for StencilMatrix --- psydac/feec/global_projectors.py | 3 ++- .../feec/tests/test_commuting_projections.py | 22 +++++++++---------- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 74b2abca3..9af6e1597 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -13,6 +13,7 @@ from psydac.fem.vector import VectorFemSpace from psydac.ddm.cart import DomainDecomposition, CartDecomposition +from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from abc import ABCMeta, abstractmethod @@ -147,7 +148,7 @@ def __init__(self, space, nquads = None): 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) + M = StencilMatrix(V_cart, V_cart, backend=PSYDAC_BACKEND_GPYCCEL) if cell == 'I': # interpolation case diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 4ad1f9591..2251799d0 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -351,11 +351,11 @@ def test_2d_commuting_pro_1(Nel, Nq, p, bc): # solve and compare u0vec = u0.coeffs u0vec_imat = I0inv.solve(P0._rhs) - assert np.allclose(u0vec.toarray(), u0vec_imat.toarray()) + 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()) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @@ -427,11 +427,11 @@ def test_2d_commuting_pro_2(Nel, Nq, p, bc): # solve and compare u0vec = u0.coeffs u0vec_imat = I0inv.solve(P0._rhs) - assert np.allclose(u0vec.toarray(), u0vec_imat.toarray()) + 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()) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [8]) @@ -664,27 +664,27 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc): # solve and compare u0vec = u0.coeffs u0vec_imat = I0inv.solve(P0._rhs) - assert np.allclose(u0vec.toarray(), u0vec_imat.toarray()) + 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()) + assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5) #============================================================================== if __name__ == '__main__': Nel = 12 - Nq = 8 + Nq = 5 p = 3 bc = False - # test_1d_commuting_pro_1(Nel, Nq, p, bc) - # test_2d_commuting_pro_1(Nel, Nq, p, bc) - # test_2d_commuting_pro_2(Nel, Nq, p, bc) + test_1d_commuting_pro_1(20, Nq, p, bc) + test_2d_commuting_pro_1(Nel, Nq, p, bc) + test_2d_commuting_pro_2(Nel, Nq, p, bc) # test_2d_commuting_pro_3(Nel, Nq, p, bc) # test_2d_commuting_pro_4(Nel, Nq, p, bc) # test_3d_commuting_pro_1(Nel, Nq, p, bc) - test_3d_commuting_pro_2(Nel, Nq, p, bc) + # test_3d_commuting_pro_2(Nel, Nq, p, bc) # test_3d_commuting_pro_3(Nel, Nq, p, bc) # test_2d_commuting_pro_4(Nel, Nq, p, bc) From 803384aa0946ce7173eedc97af52b3c419010033 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 13 Oct 2023 09:28:06 +0200 Subject: [PATCH 4/6] Remove 3d unit tests (take too long), remove backend option (only 1d matrices) --- psydac/feec/global_projectors.py | 3 +- .../feec/tests/test_commuting_projections.py | 101 +++++++++--------- 2 files changed, 53 insertions(+), 51 deletions(-) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 9af6e1597..74b2abca3 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -13,7 +13,6 @@ from psydac.fem.vector import VectorFemSpace from psydac.ddm.cart import DomainDecomposition, CartDecomposition -from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL from abc import ABCMeta, abstractmethod @@ -148,7 +147,7 @@ def __init__(self, space, nquads = None): 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, backend=PSYDAC_BACKEND_GPYCCEL) + M = StencilMatrix(V_cart, V_cart) if cell == 'I': # interpolation case diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 2251799d0..26a468c9f 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -74,27 +74,28 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc): 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, verbose=True) - I1inv = Inverse(imat_kronecker_P1, verbose=True) + # imat_kronecker_P0 = P0.imat_kronecker + # imat_kronecker_P1 = P1.imat_kronecker + # I0inv = Inverse(imat_kronecker_P0, verbose=True) + # I1inv = Inverse(imat_kronecker_P1, verbose=True) - # build the rhs - P0.func(fun1) - P1.func(D1fun1, D2fun1, D3fun1) + # # 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) + # # 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) + # 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]) @@ -171,27 +172,28 @@ def test_3d_commuting_pro_2(Nel, Nq, p, bc): 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, verbose=True) - I2inv = Inverse(imat_kronecker_P2, verbose=True) + # imat_kronecker_P1 = P1.imat_kronecker + # imat_kronecker_P2 = P2.imat_kronecker + # I1inv = Inverse(imat_kronecker_P1, verbose=True) + # I2inv = Inverse(imat_kronecker_P2, verbose=True) - # build the rhs - P1.func(fun1, fun2, fun3) - P2.func(cf1, cf2, cf3) + # # 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) + # # 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) + # 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]) @@ -259,27 +261,28 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc): 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, verbose=True) - I3inv = Inverse(imat_kronecker_P3, verbose=True) + # imat_kronecker_P2 = P2.imat_kronecker + # imat_kronecker_P3 = P3.imat_kronecker + # I2inv = Inverse(imat_kronecker_P2, verbose=True) + # I3inv = Inverse(imat_kronecker_P3, verbose=True) - # build the rhs - P2.func(fun1, fun2, fun3) - P3.func(difun) + # # 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) + # # 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) + # u3vec = u3.coeffs + # u3vec_imat = I3inv.solve(P3._rhs) + # assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-5) @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @@ -678,14 +681,14 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc): p = 3 bc = False - test_1d_commuting_pro_1(20, Nq, p, bc) + test_1d_commuting_pro_1(Nel, Nq, p, bc) test_2d_commuting_pro_1(Nel, Nq, p, bc) test_2d_commuting_pro_2(Nel, Nq, p, bc) - # test_2d_commuting_pro_3(Nel, Nq, p, bc) - # test_2d_commuting_pro_4(Nel, Nq, p, bc) - # test_3d_commuting_pro_1(Nel, Nq, p, bc) - # test_3d_commuting_pro_2(Nel, Nq, p, bc) - # test_3d_commuting_pro_3(Nel, Nq, p, bc) - # test_2d_commuting_pro_4(Nel, Nq, p, bc) + test_2d_commuting_pro_3(Nel, Nq, p, bc) + test_2d_commuting_pro_4(Nel, Nq, p, bc) + test_3d_commuting_pro_1(Nel, Nq, p, bc) + test_3d_commuting_pro_2(Nel, Nq, p, bc) + test_3d_commuting_pro_3(Nel, Nq, p, bc) + test_2d_commuting_pro_4(Nel, Nq, p, bc) From 8f7e477d572b0f1ac1a8aa9523c6af10bb3bdd7f Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Fri, 3 Nov 2023 10:53:25 +0100 Subject: [PATCH 5/6] Address comments --- .../feec/tests/test_commuting_projections.py | 35 +++++++++---------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 26a468c9f..e7322b425 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -8,7 +8,7 @@ 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 GMRES as Inverse +from psydac.linalg.solvers import inverse from mpi4py import MPI import numpy as np @@ -81,8 +81,8 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc): # build the solver from the LinearOperator # imat_kronecker_P0 = P0.imat_kronecker # imat_kronecker_P1 = P1.imat_kronecker - # I0inv = Inverse(imat_kronecker_P0, verbose=True) - # I1inv = Inverse(imat_kronecker_P1, verbose=True) + # I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + # I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) # # build the rhs # P0.func(fun1) @@ -179,8 +179,8 @@ def test_3d_commuting_pro_2(Nel, Nq, p, bc): # build the solver from the LinearOperator # imat_kronecker_P1 = P1.imat_kronecker # imat_kronecker_P2 = P2.imat_kronecker - # I1inv = Inverse(imat_kronecker_P1, verbose=True) - # I2inv = Inverse(imat_kronecker_P2, verbose=True) + # I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + # I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) # # build the rhs # P1.func(fun1, fun2, fun3) @@ -268,8 +268,8 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc): # build the solver from the LinearOperator # imat_kronecker_P2 = P2.imat_kronecker # imat_kronecker_P3 = P3.imat_kronecker - # I2inv = Inverse(imat_kronecker_P2, verbose=True) - # I3inv = Inverse(imat_kronecker_P3, verbose=True) + # I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) + # I3inv = inverse(imat_kronecker_P3, 'gmres', verbose=True) # # build the rhs # P2.func(fun1, fun2, fun3) @@ -344,8 +344,8 @@ def test_2d_commuting_pro_1(Nel, Nq, p, bc): # build the solver from the LinearOperator imat_kronecker_P0 = P0.imat_kronecker imat_kronecker_P1 = P1.imat_kronecker - I0inv = Inverse(imat_kronecker_P0, verbose=True) - I1inv = Inverse(imat_kronecker_P1, verbose=True) + I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) # build the rhs P0.func(fun1) @@ -420,8 +420,8 @@ def test_2d_commuting_pro_2(Nel, Nq, p, bc): # build the solver from the LinearOperator imat_kronecker_P0 = P0.imat_kronecker imat_kronecker_P1 = P1.imat_kronecker - I0inv = Inverse(imat_kronecker_P0, verbose=True) - I1inv = Inverse(imat_kronecker_P1, verbose=True) + I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) # build the rhs P0.func(fun1) @@ -503,8 +503,8 @@ def test_2d_commuting_pro_3(Nel, Nq, p, bc): # build the solver from the LinearOperator imat_kronecker_P2 = P2.imat_kronecker imat_kronecker_P3 = P3.imat_kronecker - I2inv = Inverse(imat_kronecker_P2, verbose=True) - I3inv = Inverse(imat_kronecker_P3, verbose=True) + I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) + I3inv = inverse(imat_kronecker_P3, 'gmres', verbose=True) # build the rhs P2.func(fun1, fun2) @@ -586,8 +586,8 @@ def test_2d_commuting_pro_4(Nel, Nq, p, bc): # build the solver from the LinearOperator imat_kronecker_P1 = P1.imat_kronecker imat_kronecker_P2 = P2.imat_kronecker - I1inv = Inverse(imat_kronecker_P1, verbose=True) - I2inv = Inverse(imat_kronecker_P2, verbose=True) + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) + I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True) # build the rhs P1.func(fun1, fun2) @@ -657,8 +657,8 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc): # build the solver from the LinearOperator imat_kronecker_P0 = P0.imat_kronecker imat_kronecker_P1 = P1.imat_kronecker - I0inv = Inverse(imat_kronecker_P0, verbose=True) - I1inv = Inverse(imat_kronecker_P1, verbose=True) + I0inv = inverse(imat_kronecker_P0, 'gmres', verbose=True) + I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) # build the rhs P0.func(fun1) @@ -689,6 +689,5 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc): test_3d_commuting_pro_1(Nel, Nq, p, bc) test_3d_commuting_pro_2(Nel, Nq, p, bc) test_3d_commuting_pro_3(Nel, Nq, p, bc) - test_2d_commuting_pro_4(Nel, Nq, p, bc) From 489130fd90610054d7dc004f6b8686568146fea7 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 14 Dec 2023 15:48:03 +0100 Subject: [PATCH 6/6] Run the 2d commuting projectors tests in parallel. --- psydac/feec/tests/test_commuting_projections.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index e7322b425..00e0e4bf4 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -284,6 +284,7 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc): # 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]) @@ -360,6 +361,7 @@ def test_2d_commuting_pro_1(Nel, Nq, p, bc): 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]) @@ -436,6 +438,7 @@ def test_2d_commuting_pro_2(Nel, Nq, p, bc): 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]) @@ -519,6 +522,7 @@ def test_2d_commuting_pro_3(Nel, Nq, p, bc): 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])