From ba29e906917670822aa2aa5264b6ce58ea692933 Mon Sep 17 00:00:00 2001 From: Stefan Possanner Date: Thu, 12 Oct 2023 15:30:38 +0200 Subject: [PATCH 01/24] 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 02/24] 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 03/24] 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 04/24] 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 05/24] 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 06/24] 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]) From 9516a403f3159b49c126f8b2c619415bee9a80aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 26 Aug 2024 15:54:17 +0200 Subject: [PATCH 07/24] Print pytest's detailed summary report --- .github/workflows/continuous-integration.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 7724f4be4..690e5c107 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -173,28 +173,28 @@ jobs: run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python -m pytest -n auto --pyargs psydac -m "not parallel and not petsc" + python -m pytest -n auto --pyargs psydac -m "not parallel and not petsc" -ra - name: Run MPI tests with Pytest working-directory: ./pytest run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and not petsc" + python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and not petsc" -ra - name: Run single-process PETSc tests with Pytest working-directory: ./pytest run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python -m pytest -n auto --pyargs psydac -m "not parallel and petsc" + python -m pytest -n auto --pyargs psydac -m "not parallel and petsc" -ra - name: Run MPI PETSc tests with Pytest - working-directory: ./pytest + working-directory: ./ run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and petsc" + python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and petsc" -ra - name: Remove test directory if: ${{ always() }} From a6ce393e599e50c50c3df56b42b344b5346b9331 Mon Sep 17 00:00:00 2001 From: Stefan Possanner <86720346+spossann@users.noreply.github.com> Date: Fri, 13 Sep 2024 14:48:38 +0200 Subject: [PATCH 08/24] Update psydac/feec/tests/test_commuting_projections.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- .../feec/tests/test_commuting_projections.py | 30 +++++++------------ 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index b2f89906a..099fad49a 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -76,28 +76,20 @@ 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) + Id_0 = IdentityLinearOperator(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 = sqrt(e0 @ e0) + assert norm2_e0 < 1e-12 + + Id_1 = IdentityLinearOperator(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 = sqrt(e1 @ e1) + assert norm2_e1 < 1e-12 #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) From 214079d70194a8a84a04282c46c65e4a5bc7320f Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Wed, 18 Sep 2024 16:50:35 +0200 Subject: [PATCH 09/24] Fixes in `test_3d_commuting_pro_1` (#438) Change `IdentityLinearOperator` to `IdentityOperator`, and `e0 @ e0` to `e0.dot(e0)`. Add missing import statement. --- psydac/feec/tests/test_commuting_projections.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 099fad49a..85054461d 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -9,6 +9,7 @@ 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 @@ -79,16 +80,16 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc, m): #-------------------------- # check BlockLinearOperator #-------------------------- - Id_0 = IdentityLinearOperator(H1.vector_space) + 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 = sqrt(e0 @ e0) + norm2_e0 = np.sqrt(e0.dot(e0)) assert norm2_e0 < 1e-12 - Id_1 = IdentityLinearOperator(Hcurl.vector_space) + 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 = sqrt(e1 @ e1) + norm2_e1 = np.sqrt(e1.dot(e1)) assert norm2_e1 < 1e-12 #============================================================================== From 245a45e06963a86c0ab31eeeaf13a84e5ce9fabc Mon Sep 17 00:00:00 2001 From: Max Date: Thu, 19 Sep 2024 14:06:51 +0200 Subject: [PATCH 10/24] Updated GlobalProjector so the m>1 is handled correctly --- psydac/feec/global_projectors.py | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 1ab23ec91..68d4749b0 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -139,13 +139,14 @@ def __init__(self, space, nquads = None): 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], [1]) + cart_decomp = CartDecomposition(domain_decomp, [n], [[s]], [[e]], [p], [m]) V_cart = StencilVectorSpace(cart_decomp) M = StencilMatrix(V_cart, V_cart) @@ -169,11 +170,12 @@ def __init__(self, space, nquads = None): 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] + 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]) - + if m == 1: + 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': @@ -202,10 +204,12 @@ def __init__(self, space, nquads = None): 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] + 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.hmat[s:e + 1]) + if m == 1: + assert np.allclose(M.toarray()[s:e + 1], V.imat[s:e + 1]) + # TODO Fix toarray() for multiplicity m > 1 matrixcells += [M.copy()] From e66f86371cb603e717acf6cd70e7d1835f5e8c9e Mon Sep 17 00:00:00 2001 From: Max Lindqvist Date: Fri, 20 Sep 2024 11:55:32 +0200 Subject: [PATCH 11/24] Fixed bug, replaced imat with hmat --- psydac/feec/global_projectors.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 68d4749b0..91f0c2569 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -204,11 +204,11 @@ def __init__(self, space, nquads = None): 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] + 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 if m == 1: - assert np.allclose(M.toarray()[s:e + 1], V.imat[s:e + 1]) + assert np.allclose(M.toarray()[s:e + 1], V.hmat[s:e + 1]) # TODO Fix toarray() for multiplicity m > 1 matrixcells += [M.copy()] From cd8e480fc773bb601a5de2d3583a97084db3a430 Mon Sep 17 00:00:00 2001 From: Max Date: Mon, 23 Sep 2024 15:45:15 +0200 Subject: [PATCH 12/24] Fixed indexing in stencil2coo_1d_C adn stencil2coo_1d_F for shifts > 1, this fixes the problem with toarray() which did previously not work for StencilMatrix. Updated indexing in KroneckerStencilMatrix to include shift factors in index calculations. Added back the assert statement in GlobalProjector since the toarray() now works. The kernels still don't work properly for different shifts in the domain and the co-domain, but I believe equivalent changes can be done in the remaining kernels. --- psydac/feec/global_projectors.py | 7 ++----- psydac/linalg/kernels/stencil2coo_kernels.py | 12 +++++++++--- psydac/linalg/kron.py | 8 ++++---- 3 files changed, 15 insertions(+), 12 deletions(-) diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 91f0c2569..509995406 100644 --- a/psydac/feec/global_projectors.py +++ b/psydac/feec/global_projectors.py @@ -173,8 +173,7 @@ def __init__(self, space, nquads = None): 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 - if m == 1: - assert np.allclose(M.toarray()[s:e + 1], V.imat[s:e + 1]) + assert np.allclose(M.toarray()[s:e + 1], V.imat[s:e + 1]) # TODO Fix toarray() for multiplicity m > 1 matrixcells += [M.copy()] @@ -207,9 +206,7 @@ def __init__(self, space, nquads = None): 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 - if m == 1: - assert np.allclose(M.toarray()[s:e + 1], V.hmat[s:e + 1]) - # TODO Fix toarray() for multiplicity m > 1 + assert np.allclose(M.toarray()[s:e + 1], V.hmat[s:e + 1]) matrixcells += [M.copy()] diff --git a/psydac/linalg/kernels/stencil2coo_kernels.py b/psydac/linalg/kernels/stencil2coo_kernels.py index 55f4fc874..3ca953997 100644 --- a/psydac/linalg/kernels/stencil2coo_kernels.py +++ b/psydac/linalg/kernels/stencil2coo_kernels.py @@ -3,9 +3,9 @@ #!!!!!!!!!!!!! #TODO avoid using The expensive modulo operator % in the non periodic case to make the methods faster +#TODO The kernels don't work properly for different shifts in the domain and the codomain #!!!!!!!!!!!!! - #__all__ = ['stencil2coo_1d_C','stencil2coo_1d_F','stencil2coo_2d_C','stencil2coo_2d_F', 'stencil2coo_3d_C', 'stencil2coo_3d_F'] #======================================================================================================== @@ -19,7 +19,10 @@ def stencil2coo_1d_C(A:'T[:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:]', for j1 in range(ncl1): value = A[i1+pp1,j1] if abs(value) == 0.0:continue - J = ((I//cm1)*dm1+j1-dp1)%nc1 + if cm1 == dm1: + J = (I + j1 - dp1)%nc1 + else: + J = ((I//cm1)*dm1+j1-dp1)%nc1 rows[nnz] = I cols[nnz] = J data[nnz] = value @@ -38,7 +41,10 @@ def stencil2coo_1d_F(A:'T[:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:]', for j1 in range(ncl1): value = A[i1+pp1,j1] if abs(value) == 0.0:continue - J = ((I//cm1)*dm1+j1-dp1)%nc1 + if cm1 == dm1: + J = (I + j1 - dp1)%nc1 + else: + J = ((I//cm1)*dm1+j1-dp1)%nc1 rows[nnz] = I cols[nnz] = J data[nnz] = value 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 From 28e0f89b2b72a07e461cf842b482dcfee5a6be63 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 6 Sep 2024 15:41:45 +0200 Subject: [PATCH 13/24] Improve and test Pyccel compiler flags (#428) - Add compiler flags for GFortran >= 14 on Apple silicon chips - Add unit test in `test_epyccel_flags.py` to check compiler flags passed to `epyccel` - Update GitHub actions used in documentation workflow --- .github/workflows/continuous-integration.yml | 4 +-- .github/workflows/documentation.yml | 10 +++--- psydac/api/settings.py | 19 ++++++++--- psydac/api/tests/test_epyccel_flags.py | 35 ++++++++++++++++++++ pytest.ini | 1 + 5 files changed, 58 insertions(+), 11 deletions(-) create mode 100644 psydac/api/tests/test_epyccel_flags.py diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 690e5c107..557bd570e 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -159,9 +159,9 @@ jobs: python -m pip install . python -m pip freeze - - name: Test pyccel optimization flags + - name: Test Pyccel optimization flags run: | - pyccel --verbose ./psydac/linalg/kernels/axpy_kernels.py + pytest --pyargs psydac -m pyccel --capture=no - name: Initialize test directory run: | diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 88e8065d5..be24f89c9 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -18,9 +18,9 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN}} steps: - name: Checkout - uses: actions/checkout@v3 + uses: actions/checkout@v4 - name: Set up Python - uses: actions/setup-python@v4 + uses: actions/setup-python@v5 with: python-version: 3.9 - name: Install non-Python dependencies on Ubuntu @@ -37,9 +37,9 @@ jobs: make -C docs html python docs/update_links.py - name: Setup Pages - uses: actions/configure-pages@v3 + uses: actions/configure-pages@v5 - name: Upload artifact - uses: actions/upload-pages-artifact@v1 + uses: actions/upload-pages-artifact@v3 with: path: 'docs/build/html' @@ -53,4 +53,4 @@ jobs: steps: - name: Deploy to GitHub Pages id: deployment - uses: actions/deploy-pages@v2 + uses: actions/deploy-pages@v4 diff --git a/psydac/api/settings.py b/psydac/api/settings.py index 74a48c645..3a3830107 100644 --- a/psydac/api/settings.py +++ b/psydac/api/settings.py @@ -1,5 +1,8 @@ # coding: utf-8 +import subprocess # nosec B404 import platform +import re +from packaging.version import Version __all__ = ('PSYDAC_DEFAULT_FOLDER', 'PSYDAC_BACKENDS') @@ -40,18 +43,26 @@ 'openmp' : False} # ... +# Get gfortran version +gfortran_version_output = subprocess.check_output(['gfortran', '--version']).decode('utf-8') # nosec B603, B607 +gfortran_version_string = re.search("(\d+\.\d+\.\d+)", gfortran_version_output).group() +gfortran_version = Version(gfortran_version_string) + # Platform-dependent flags -if platform.system() == "Darwin" and platform.machine() == 'arm64': +if platform.system() == "Darwin" and platform.machine() == 'arm64' and gfortran_version >= Version("14"): + # Apple silicon requires architecture-specific flags (see https://github.com/pyccel/psydac/pull/411) - import subprocess # nosec B404 + # which are only available on GCC version >= 14 cpu_brand = subprocess.check_output(['sysctl','-n','machdep.cpu.brand_string']).decode('utf-8') # nosec B603, B607 - if "Apple M1" in cpu_brand: - PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m1' + if "Apple M1" in cpu_brand: PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m1' + elif "Apple M2" in cpu_brand: PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m2' + elif "Apple M3" in cpu_brand: PSYDAC_BACKEND_GPYCCEL['flags'] += ' -mcpu=apple-m3' else: # TODO: Support later Apple CPU models. Perhaps the CPU naming scheme could be easily guessed # based on the output of 'sysctl -n machdep.cpu.brand_string', but I wouldn't rely on this # guess unless it has been manually verified. Loud errors are better than silent failures! raise SystemError(f"Unsupported Apple CPU '{cpu_brand}'.") + else: # Default architecture flags PSYDAC_BACKEND_GPYCCEL['flags'] += ' -march=native -mtune=native' diff --git a/psydac/api/tests/test_epyccel_flags.py b/psydac/api/tests/test_epyccel_flags.py new file mode 100644 index 000000000..cdbcb0906 --- /dev/null +++ b/psydac/api/tests/test_epyccel_flags.py @@ -0,0 +1,35 @@ +import pytest + + +@pytest.mark.pyccel +def test_epyccel_flags(): + + from pyccel.epyccel import epyccel + from psydac.api.settings import PSYDAC_BACKEND_GPYCCEL as backend + + kwargs = {'language' : 'fortran', + 'compiler' : backend['compiler'], + 'fflags' : backend['flags'], + 'accelerators' : ['openmp'] if backend['openmp'] else [], + 'verbose' : True, + } + + # Function to be Pyccel-ized + def f(x : float): + return 3 * x + + # Pyccel magic + # ------------ + # Fortran code is generated and then compiled with the selected compiler. + # The compiled function is callable from Python through the C Python API. + # The necessary C wrapper functions are also generated by Pyccel. + fast_f = epyccel(f, **kwargs) + + # Check output of pyccelized function + assert fast_f(3.5) == f(3.5) + + +# Interactive usage +if __name__ == '__main__': + test_epyccel_flags() + print('PASSED') diff --git a/pytest.ini b/pytest.ini index c18eb444d..91ef13cfa 100644 --- a/pytest.ini +++ b/pytest.ini @@ -9,6 +9,7 @@ markers = serial: single-process test, parallel: test to be run using 'mpiexec', petsc: test requiring a working PETSc installation with petsc4py Python bindings + pyccel: test for checking Pyccel setup on machine python_files = test_*.py python_classes = From dd26dd56d47480a7a2bb94ff2b440e8aea19d789 Mon Sep 17 00:00:00 2001 From: Frederik Schnack Date: Wed, 18 Sep 2024 16:42:46 +0200 Subject: [PATCH 14/24] Update calls to `numpy.reshape` to avoid deprecation warnings (#433) NumPy 2.1 deprecates the `newshape` argument of the `reshape` function. Instead, one should use the positional-only argument `shape`. Every call in Psydac has been updated accordingly. This closes #432. --- psydac/fem/tensor.py | 2 +- psydac/mapping/discrete.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py index 8832ad1a8..e0367e43d 100644 --- a/psydac/fem/tensor.py +++ b/psydac/fem/tensor.py @@ -456,7 +456,7 @@ def eval_fields(self, grid, *fields, weights=None, npts_per_cell=None, overlap=0 npts_per_cell = (npts_per_cell,) * self.ldim for i in range(self.ldim): ncells_i = len(self.breaks[i]) - 1 - grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i])) + grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i])) out_fields = self.eval_fields_regular_tensor_grid(grid, *fields, weights=weights, overlap=overlap) # return a list return [np.ascontiguousarray(out_fields[..., i]) for i in range(len(fields))] diff --git a/psydac/mapping/discrete.py b/psydac/mapping/discrete.py index 636b46b0f..af4c54ba3 100644 --- a/psydac/mapping/discrete.py +++ b/psydac/mapping/discrete.py @@ -243,7 +243,7 @@ def jac_mat_grid(self, grid, npts_per_cell=None, overlap=0): npts_per_cell = (npts_per_cell,) * self.ldim for i in range(self.ldim): ncells_i = len(self.space.breaks[i]) - 1 - grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i])) + grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i])) jac_mats = self.jac_mat_regular_tensor_grid(grid, overlap=overlap) return jac_mats @@ -413,7 +413,7 @@ def inv_jac_mat_grid(self, grid, npts_per_cell=None, overlap=0): npts_per_cell = (npts_per_cell,) * self.ldim for i in range(self.ldim): ncells_i = len(self.space.breaks[i]) - 1 - grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i])) + grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i])) inv_jac_mats = self.inv_jac_mat_regular_tensor_grid(grid, overlap=overlap) return inv_jac_mats @@ -584,7 +584,7 @@ def jac_det_grid(self, grid, npts_per_cell=None, overlap=0): npts_per_cell = (npts_per_cell,) * self.ldim for i in range(self.ldim): ncells_i = len(self.space.breaks[i]) - 1 - grid[i] = np.reshape(grid[i], newshape=(ncells_i, npts_per_cell[i])) + grid[i] = np.reshape(grid[i], (ncells_i, npts_per_cell[i])) jac_dets = self.jac_det_regular_tensor_grid(grid, overlap=overlap) return jac_dets From a3df57d4d75ba4df165b037bd14e5fe609b09f8b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 19 Sep 2024 10:46:18 +0200 Subject: [PATCH 15/24] Circumvent `apt-get` download error (#440) Update `apt` index files with `apt-get update` to avoid download error in CI workflow. --- .github/workflows/continuous-integration.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 557bd570e..8b6dc33cf 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -58,6 +58,7 @@ jobs: - name: Set default MPI and HDF5 C compilers on Ubuntu if: matrix.os == 'ubuntu-latest' run: | + sudo apt-get update sudo apt-get install --reinstall openmpi-bin libhdf5-openmpi-dev - name: Install non-Python dependencies on macOS From 188f2d228e7103f9d4baa271174778aa2560b956 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 19 Sep 2024 12:07:41 +0200 Subject: [PATCH 16/24] Mention C support in near future in README (#441) Clarify in the `README.md` file that, although the C backend may be selected for accelerating the kernel files with Pyccel, this is not fully working yet. Hence the Fortran backend (which is the default) is the only one available. A future version of Pyccel will certainly provide a C backend as capable as the Fortran one. See issue #431. Co-authored-by: Martin Campos Pinto --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index c35348982..8948a5efd 100644 --- a/README.md +++ b/README.md @@ -186,7 +186,7 @@ python3 /mpi_tester.py --pyargs psydac -m "parallel and petsc" Many of Psydac's low-level Python functions can be translated to a compiled language using the [Pyccel](https://github.com/pyccel/pyccel) transpiler. Currently, all of those functions are collected in modules which follow the name pattern `[module]_kernels.py`. -The classical installation translates all kernel files to Fortran without user intervention. This does not happen in the case of an editable install, but the command `psydac-accelerate` is made available to the user instead. This command applies Pyccel to all the kernel files in the source directory, and the C language may be selected instead of Fortran (which is the default). +The classical installation translates all kernel files to Fortran without user intervention. This does not happen in the case of an editable install, but the command `psydac-accelerate` is made available to the user instead. This command applies Pyccel to all the kernel files in the source directory. The default language is currently Fortran, C should also be supported in a near future. - **Only in development mode**: ```bash From c617b15979454f6c11514bd6b05709b7d9c224de Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Wed, 25 Sep 2024 14:46:06 +0200 Subject: [PATCH 17/24] Remove `is_block` property from `VectorFemSpace` (#439) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The obsolete property `is_block` was never used in Psydac nor in Struphy, hence we remove it. --------- Co-authored-by: Yaman Güçlü --- psydac/fem/vector.py | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/psydac/fem/vector.py b/psydac/fem/vector.py index bca983331..84a29db10 100644 --- a/psydac/fem/vector.py +++ b/psydac/fem/vector.py @@ -333,21 +333,6 @@ def ncells(self): def spaces( self ): return self._spaces - @property - def is_block(self): - """Returns True if all components are identical spaces.""" - # TODO - improve this tests. for the moment, we only check the degree, - # - shall we check the bc too? - - degree = [V.degree for V in self.spaces] - if self.pdim == 1: - return len(np.unique(degree)) == 1 - else: - ns = np.asarray(degree[0]) - for ms in degree[1:]: - if not( np.allclose(ns, np.asarray(ms)) ): return False - return True - # ... def get_refined_space(self, ncells): return self._refined_space[tuple(ncells)] From 6fc53bbcd5cef3d26cc135b42e18dd5067ac14d5 Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Fri, 27 Sep 2024 16:24:41 +0200 Subject: [PATCH 18/24] Update `stencil2coo` kernels for shifts > 1 (#442) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fix the indexing in all kernels in `stencil2coo_kernels.py` so that the method `toarray` works for `StencilMatrix` objects in 1D/2D/3D with `shift` larger than 1. The unit tests for `toarray` in `test_stencil_matrix.py` have been parametrized to run with `shift` values of 1 and 2. --------- Co-authored-by: Yaman Güçlü --- psydac/linalg/kernels/stencil2coo_kernels.py | 30 ++++++++------------ psydac/linalg/tests/test_stencil_matrix.py | 8 +++--- 2 files changed, 16 insertions(+), 22 deletions(-) diff --git a/psydac/linalg/kernels/stencil2coo_kernels.py b/psydac/linalg/kernels/stencil2coo_kernels.py index 3ca953997..464fcab8e 100644 --- a/psydac/linalg/kernels/stencil2coo_kernels.py +++ b/psydac/linalg/kernels/stencil2coo_kernels.py @@ -19,10 +19,7 @@ def stencil2coo_1d_C(A:'T[:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:]', for j1 in range(ncl1): value = A[i1+pp1,j1] if abs(value) == 0.0:continue - if cm1 == dm1: - J = (I + j1 - dp1)%nc1 - else: - J = ((I//cm1)*dm1+j1-dp1)%nc1 + J = ((I*dm1//cm1)+j1-dp1)%nc1 rows[nnz] = I cols[nnz] = J data[nnz] = value @@ -41,10 +38,7 @@ def stencil2coo_1d_F(A:'T[:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[:]', for j1 in range(ncl1): value = A[i1+pp1,j1] if abs(value) == 0.0:continue - if cm1 == dm1: - J = (I + j1 - dp1)%nc1 - else: - J = ((I//cm1)*dm1+j1-dp1)%nc1 + J = ((I*dm1//cm1)+j1-dp1)%nc1 rows[nnz] = I cols[nnz] = J data[nnz] = value @@ -72,8 +66,8 @@ def stencil2coo_2d_C(A:'T[:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[: for j2 in range(ncl2): value = A[i1+pp1,i2+pp2,j1,j2] if abs(value) == 0.0:continue - jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1 - jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2 + jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1 + jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2 J = jj1*nc2 + jj2 @@ -103,8 +97,8 @@ def stencil2coo_2d_F(A:'T[:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int64[: for j2 in range(ncl2): value = A[i1+pp1,i2+pp2,j1,j2] if abs(value) == 0.0:continue - jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1 - jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2 + jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1 + jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2 J = jj2*nc1 + jj1 @@ -138,9 +132,9 @@ def stencil2coo_3d_C(A:'T[:,:,:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int for j3 in range(ncl3): value = A[i1+pp1,i2+pp2,i3+pp3,j1,j2,j3] if abs(value) == 0.0:continue - jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1 - jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2 - jj3 = ((ii3//cm3)*dm3+j3-dp3)%nc3 + jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1 + jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2 + jj3 = ((ii3*dm3//cm3)+j3-dp3)%nc3 J = jj1*nc2*nc3 + jj2*nc3 + jj3 @@ -176,9 +170,9 @@ def stencil2coo_3d_F(A:'T[:,:,:,:,:,:]', data:'T[:]', rows:'int64[:]', cols:'int for j3 in range(ncl3): value = A[i1+pp1,i2+pp2,i3+pp3,j1,j2,j3] if abs(value) == 0.0:continue - jj1 = ((ii1//cm1)*dm1+j1-dp1)%nc1 - jj2 = ((ii2//cm2)*dm2+j2-dp2)%nc2 - jj3 = ((ii3//cm3)*dm3+j3-dp3)%nc3 + jj1 = ((ii1*dm1//cm1)+j1-dp1)%nc1 + jj2 = ((ii2*dm2//cm2)+j2-dp2)%nc2 + jj3 = ((ii3*dm3//cm3)+j3-dp3)%nc3 J = jj3*nc1*nc2 + jj2*nc1 + jj1 diff --git a/psydac/linalg/tests/test_stencil_matrix.py b/psydac/linalg/tests/test_stencil_matrix.py index 404a0f970..cf2e602a6 100644 --- a/psydac/linalg/tests/test_stencil_matrix.py +++ b/psydac/linalg/tests/test_stencil_matrix.py @@ -433,7 +433,7 @@ def test_stencil_matrix_2d_serial_spurious_entries( dtype, p1, p2, s1, s2, P1, P @pytest.mark.parametrize('dtype', [float]) @pytest.mark.parametrize('n1', [7]) @pytest.mark.parametrize('p1', [1,2]) -@pytest.mark.parametrize('s1', [1]) +@pytest.mark.parametrize('s1', [1, 2]) @pytest.mark.parametrize('P1', [True, False]) def test_stencil_matrix_1d_serial_toarray( dtype, n1, p1, s1, P1): # Select non-zero values based on diagonal index @@ -489,8 +489,8 @@ def test_stencil_matrix_1d_serial_toarray( dtype, n1, p1, s1, P1): @pytest.mark.parametrize('n2', [8, 7]) @pytest.mark.parametrize('p1', [1, 2]) @pytest.mark.parametrize('p2', [1, 3]) -@pytest.mark.parametrize('s1', [1]) -@pytest.mark.parametrize('s2', [1]) +@pytest.mark.parametrize('s1', [1, 2]) +@pytest.mark.parametrize('s2', [1, 2]) @pytest.mark.parametrize('P1', [True, False]) @pytest.mark.parametrize('P2', [True, False]) def test_stencil_matrix_2d_serial_toarray( dtype, n1, n2, p1, p2, s1, s2, P1, P2): @@ -2191,7 +2191,7 @@ def test_stencil_matrix_2d_serial_backend_switch(dtype, n1, n2, p1, p2, s1, s2, @pytest.mark.parametrize('dtype', [float, complex]) @pytest.mark.parametrize('n1', [20, 67]) @pytest.mark.parametrize('p1', [1, 2, 3]) -@pytest.mark.parametrize('sh1', [1]) +@pytest.mark.parametrize('sh1', [1, 2]) @pytest.mark.parametrize('P1', [True, False]) @pytest.mark.parallel def test_stencil_matrix_1d_parallel_toarray(dtype, n1, p1, sh1, P1): From 2dc2297a76e9d3372819d22371d440847733a45a Mon Sep 17 00:00:00 2001 From: maxlin-ipp Date: Fri, 27 Sep 2024 17:37:55 +0200 Subject: [PATCH 19/24] Add `sqrt` option to `diagonal` method of stencil/block matrices (#434) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --------- Co-authored-by: Yaman Güçlü --- psydac/linalg/block.py | 9 +++++++-- psydac/linalg/stencil.py | 10 +++++++++- psydac/linalg/tests/test_linalg.py | 6 ++++++ 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index d54e63626..49e29b60d 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -777,13 +777,18 @@ def __sub__(self, M): #-------------------------------------- # New properties/methods #-------------------------------------- - def diagonal(self, *, inverse = False, out = None): + def diagonal(self, *, inverse = False, sqrt = False, out = None): """Get the coefficients on the main diagonal as another BlockLinearOperator object. Parameters ---------- inverse : bool If True, get the inverse of the diagonal. (Default: False). + Can be combined with sqrt to get the inverse square root. + + sqrt : bool + If True, get the square root of the diagonal. (Default: False). + Can be combined with inverse to get the inverse square root. out : BlockLinearOperator If provided, write the diagonal entries into this matrix. (Default: None). @@ -815,7 +820,7 @@ def diagonal(self, *, inverse = False, out = None): # Store the diagonal (or its inverse) into `out` for i, j in self.nonzero_block_indices: if i == j: - out[i, i] = self[i, i].diagonal(inverse = inverse, out = out[i, i]) + out[i, i] = self[i, i].diagonal(inverse = inverse, sqrt = sqrt, out = out[i, i]) return out diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 0c818af26..882d2c854 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -1372,7 +1372,7 @@ def _exchange_assembly_data_serial(self): self._data[idx_to] += self._data[idx_from] # ... - def diagonal(self, *, inverse = False, out = None): + def diagonal(self, *, inverse = False, sqrt = False, out = None): """ Get the coefficients on the main diagonal as a StencilDiagonalMatrix object. @@ -1380,6 +1380,11 @@ def diagonal(self, *, inverse = False, out = None): ---------- inverse : bool If True, get the inverse of the diagonal. (Default: False). + Can be combined with sqrt to get the inverse square root. + + sqrt : bool + If True, get the square root of the diagonal. (Default: False). + Can be combined with inverse to get the inverse square root. out : StencilDiagonalMatrix If provided, write the diagonal entries into this matrix. (Default: None). @@ -1418,6 +1423,9 @@ def diagonal(self, *, inverse = False, out = None): else: data = diag.copy() + if sqrt: + np.sqrt(data, out=data) + # If needed create a new StencilDiagonalMatrix object if out is None: out = StencilDiagonalMatrix(V, W, data) diff --git a/psydac/linalg/tests/test_linalg.py b/psydac/linalg/tests/test_linalg.py index b34e222c1..c9962c616 100644 --- a/psydac/linalg/tests/test_linalg.py +++ b/psydac/linalg/tests/test_linalg.py @@ -618,6 +618,12 @@ def test_inverse_transpose_interaction(n1, n2, p1, p2, P1=False, P2=False): ### -1,T & T,-1 --- -1,T,T --- -1,T,-1 --- T,-1,-1 --- T,-1,T (the combinations I test) ### + # Square root test + scaled_matrix = B * np.random.random() # Ensure the diagonal elements != 1 + diagonal_values = scaled_matrix.diagonal(sqrt=False).toarray() + sqrt_diagonal_values = scaled_matrix.diagonal(sqrt=True).toarray() + assert np.array_equal(sqrt_diagonal_values, np.sqrt(diagonal_values)) + tol = 1e-5 C = inverse(B, 'cg', tol=tol) P = B.diagonal(inverse=True) From 6783d3fd86a99f7882f23fdd258d8d180b88347e Mon Sep 17 00:00:00 2001 From: Max Date: Tue, 1 Oct 2024 15:35:42 +0200 Subject: [PATCH 20/24] reverted stencil2coo_kernels.py to devel version --- psydac/linalg/kernels/stencil2coo_kernels.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/kernels/stencil2coo_kernels.py b/psydac/linalg/kernels/stencil2coo_kernels.py index 464fcab8e..81c278df6 100644 --- a/psydac/linalg/kernels/stencil2coo_kernels.py +++ b/psydac/linalg/kernels/stencil2coo_kernels.py @@ -3,9 +3,9 @@ #!!!!!!!!!!!!! #TODO avoid using The expensive modulo operator % in the non periodic case to make the methods faster -#TODO The kernels don't work properly for different shifts in the domain and the codomain #!!!!!!!!!!!!! + #__all__ = ['stencil2coo_1d_C','stencil2coo_1d_F','stencil2coo_2d_C','stencil2coo_2d_F', 'stencil2coo_3d_C', 'stencil2coo_3d_F'] #======================================================================================================== From f75a3335f36758d23243ee5dfb251f1a261df632 Mon Sep 17 00:00:00 2001 From: Max Date: Wed, 9 Oct 2024 08:50:50 +0200 Subject: [PATCH 21/24] Updated the checks in test_commuting_projectors.py. --- .../feec/tests/test_commuting_projections.py | 220 ++++++++---------- 1 file changed, 91 insertions(+), 129 deletions(-) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 85054461d..637c50746 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -16,7 +16,7 @@ import pytest #============================================================================== -# Run test +# 3D tests #============================================================================== @pytest.mark.parametrize('Nel', [8, 12]) @pytest.mark.parametrize('Nq', [5]) @@ -92,7 +92,6 @@ def test_3d_commuting_pro_1(Nel, Nq, p, bc, m): 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]) @@ -173,26 +172,19 @@ def test_3d_commuting_pro_2(Nel, Nq, p, bc, m): #-------------------------- # 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) -#============================================================================== + 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]) @@ -264,25 +256,22 @@ def test_3d_commuting_pro_3(Nel, Nq, p, bc, m): #-------------------------- # 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) + 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]) @@ -343,24 +332,18 @@ def test_2d_commuting_pro_1(Nel, Nq, p, bc, m): #-------------------------- # 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) + + 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]) @@ -422,24 +405,18 @@ def test_2d_commuting_pro_2(Nel, Nq, p, bc, m): #-------------------------- # 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) + + 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]) @@ -508,24 +485,18 @@ def test_2d_commuting_pro_3(Nel, Nq, p, bc, m): #-------------------------- # 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) + + 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]) @@ -594,25 +565,22 @@ def test_2d_commuting_pro_4(Nel, Nq, p, bc, m): #-------------------------- # 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) + 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]) @@ -667,24 +635,18 @@ def test_1d_commuting_pro_1(Nel, Nq, p, bc, m): #-------------------------- # 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) + + 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__': From 549777dee71e5903709fed387ad47c38b7f01ca6 Mon Sep 17 00:00:00 2001 From: Max Date: Wed, 16 Oct 2024 11:55:17 +0200 Subject: [PATCH 22/24] Removed TODO from test_3d_commuting_pro_2 --- psydac/feec/tests/test_commuting_projections.py | 1 - 1 file changed, 1 deletion(-) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 637c50746..1baf6a9b9 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -168,7 +168,6 @@ 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 #-------------------------- From 8c2c79a73267304b2fbff8447a63b674e8366c29 Mon Sep 17 00:00:00 2001 From: Max Date: Wed, 16 Oct 2024 11:55:53 +0200 Subject: [PATCH 23/24] Removed TODO from test_3d_commuting_pro_3 --- psydac/feec/tests/test_commuting_projections.py | 1 - 1 file changed, 1 deletion(-) diff --git a/psydac/feec/tests/test_commuting_projections.py b/psydac/feec/tests/test_commuting_projections.py index 1baf6a9b9..006439b01 100644 --- a/psydac/feec/tests/test_commuting_projections.py +++ b/psydac/feec/tests/test_commuting_projections.py @@ -251,7 +251,6 @@ 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 #-------------------------- From f5930d160d8cc66f7a7a2f64c088fcd9ac80ca41 Mon Sep 17 00:00:00 2001 From: maxlin Date: Wed, 23 Oct 2024 16:08:55 +0200 Subject: [PATCH 24/24] Use the devel version of continuous-integration.yml --- .github/workflows/continuous-integration.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml index 8b6dc33cf..ae8d0ab61 100644 --- a/.github/workflows/continuous-integration.yml +++ b/.github/workflows/continuous-integration.yml @@ -174,28 +174,28 @@ jobs: run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python -m pytest -n auto --pyargs psydac -m "not parallel and not petsc" -ra + python -m pytest -n auto --pyargs psydac -m "not parallel and not petsc" - name: Run MPI tests with Pytest working-directory: ./pytest run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and not petsc" -ra + python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and not petsc" - name: Run single-process PETSc tests with Pytest working-directory: ./pytest run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python -m pytest -n auto --pyargs psydac -m "not parallel and petsc" -ra + python -m pytest -n auto --pyargs psydac -m "not parallel and petsc" - name: Run MPI PETSc tests with Pytest - working-directory: ./ + working-directory: ./pytest run: | export PSYDAC_MESH_DIR=$GITHUB_WORKSPACE/mesh export OMP_NUM_THREADS=2 - python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and petsc" -ra + python mpi_tester.py --mpirun="mpiexec -n 4 ${MPI_OPTS}" --pyargs psydac -m "parallel and petsc" - name: Remove test directory if: ${{ always() }}