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

spossann · 2023-10-12T13:42:41Z

Addresses issue #346.

The inter-/histopolation matrices are now exposed in GlobalProjector as KroneckerStencilMatrix or BlockLinearOperator (for vector-valued spaces).

The unit tests in feec/tests/test_commuting_projections.py have been enhanced: they now also test the new LinearOperator in GlobalProjector.imat_kronecker by computing its inverse using GMRES and comparing to the original BlockDiagonalSolver solution. 3D tests are commented out for the moment as they take too long (better solution)?

…matrices)

spossann · 2023-10-26T11:45:57Z

Hi @yguclu, it seems I cannot see why the latest test failed, do you have an idea? Also, did you have time to think about this PR yet?

psydac/feec/tests/test_commuting_projections.py

spossann · 2023-11-20T10:24:11Z

Python 3.11 tests passed upon relaunching test, without changes.

yguclu · 2023-12-05T10:06:26Z

@spossann It seems that the code could be simplified with these steps:

Use the function array_to_psydac in module psydac.linalg.utilities to convert the 1D matrices from Numpy arrays to a StencilMatrix format (not distributed);
Create a distributed 3D KroneckerStencilMatrix object from the serial 1D StencilMatrix objects;
If needed, call the method KroneckerStencilMatrix.tostencil() to obtain a distributed 3D StencilMatrix object.

Before making any of these changes, could you please add a parallel test to this PR?

spossann · 2023-12-15T07:30:50Z

@spossann It seems that the code could be simplified with these steps:

1. Use the function `array_to_psydac` in module `psydac.linalg.utilities` to convert the 1D matrices from Numpy arrays to a `StencilMatrix` format (not distributed);

The function array_to_psydac converts numpy arrays into Stencil/BlockVector only (not Matrix).

psydac/feec/tests/test_commuting_projections.py

… duplicated __del__ in StencilVector

Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

Change `IdentityLinearOperator` to `IdentityOperator`, and `e0 @ e0` to `e0.dot(e0)`. Add missing import statement.

…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.

- 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

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.

Update `apt` index files with `apt-get update` to avoid download error in CI workflow.

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 <campos@ann.jussieu.fr>

The obsolete property `is_block` was never used in Psydac nor in Struphy, hence we remove it. --------- Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

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ü <yaman.guclu@gmail.com>

--------- Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

…lation-matrices

yguclu · 2024-10-02T14:43:40Z

psydac/feec/tests/test_commuting_projections.py

+    # TODO: test takes too long in 3D
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+    # build the solver from the LinearOperator
+    # imat_kronecker_P1 = P1.imat_kronecker
+    # imat_kronecker_P2 = P2.imat_kronecker 
+    # I1inv = inverse(imat_kronecker_P1, 'gmres', verbose=True) 
+    # I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True)
+
+    # # build the rhs
+    # P1.func(fun1, fun2, fun3)
+    # P2.func(cf1, cf2, cf3)
+
+    # # solve and compare
+    # u1vec = u1.coeffs
+    # u1vec_imat = I1inv.solve(P1._rhs)
+    # assert np.allclose(u1vec.toarray(), u1vec_imat.toarray(), atol=1e-5)
+
+    # u2vec = u2.coeffs
+    # u2vec_imat = I2inv.solve(P2._rhs)
+    # assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5)
+


Why not taking the same approach I suggested in function test_3d_commuting_pro_1?

yguclu · 2024-10-02T14:44:04Z

psydac/feec/tests/test_commuting_projections.py

+    # TODO: test takes too long in 3D
+    #--------------------------
+    # check BlockLinearOperator
+    #--------------------------
+    # build the solver from the LinearOperator
+    # imat_kronecker_P2 = P2.imat_kronecker
+    # imat_kronecker_P3 = P3.imat_kronecker 
+    # I2inv = inverse(imat_kronecker_P2, 'gmres', verbose=True)
+    # I3inv = inverse(imat_kronecker_P3, 'gmres', verbose=True) 
+
+    # # build the rhs
+    # P2.func(fun1, fun2, fun3)
+    # P3.func(difun)
+
+    # # solve and compare
+    # u2vec = u2.coeffs
+    # u2vec_imat = I2inv.solve(P2._rhs)
+    # assert np.allclose(u2vec.toarray(), u2vec_imat.toarray(), atol=1e-5)
+
+    # u3vec = u3.coeffs
+    # u3vec_imat = I3inv.solve(P3._rhs)
+    # assert np.allclose(u3vec.toarray(), u3vec_imat.toarray(), atol=1e-5)


Why not taking the same approach I suggested in function test_3d_commuting_pro_1?

yguclu · 2024-10-02T14:45:01Z

psydac/feec/tests/test_commuting_projections.py

+    #--------------------------
+    # 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)


Same comment as above

yguclu · 2024-10-02T14:45:17Z

psydac/feec/tests/test_commuting_projections.py

+    #--------------------------
+    # 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)


Same comment as above

yguclu · 2024-10-02T14:45:34Z

psydac/feec/tests/test_commuting_projections.py

+    #--------------------------
+    # 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)


Same comment as above

yguclu · 2024-10-02T14:45:52Z

psydac/feec/tests/test_commuting_projections.py

+    #--------------------------
+    # 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)


Same comment as above

yguclu · 2024-10-02T14:46:05Z

psydac/feec/tests/test_commuting_projections.py

+    #--------------------------
+    # 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)


Same comment as above

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

ba29e90

spossann requested a review from yguclu October 12, 2023 13:42

camillo81 added 5 commits October 12, 2023 19:47

Set atol=1e-5 in tests

60592c9

Merge branch 'devel' into expose-inter-histopolation-matrices

47d68c6

Add forgotten atol, set pyccel backend for StencilMatrix

b59195a

Remove 3d unit tests (take too long), remove backend option (only 1d …

803384a

…matrices)

Merge branch 'devel' into expose-inter-histopolation-matrices

5fd7f65

yguclu requested changes Nov 2, 2023

View reviewed changes

psydac/feec/tests/test_commuting_projections.py Outdated Show resolved Hide resolved

psydac/feec/tests/test_commuting_projections.py Outdated Show resolved Hide resolved

Address comments

8f7e477

Merge branch 'devel' into expose-inter-histopolation-matrices

c290867

Run the 2d commuting projectors tests in parallel.

489130f

campospinto mentioned this pull request Jul 24, 2024

Expose inter-/histopolation matrices in GlobalProjector, add unit tests campospinto/psydac_dev#13

Open

Merge branch 'devel' into expose-inter-histopolation-matrices

39af2a0

yguclu reviewed Aug 22, 2024

View reviewed changes

psydac/feec/tests/test_commuting_projections.py Outdated Show resolved Hide resolved

Print pytest's detailed summary report

9516a40

spossann referenced this pull request Sep 12, 2024

Add inter-/histopolation matrices as properties of global projectors

5b3cb50

maxlin-ipp referenced this pull request Sep 12, 2024

Replace BlockMatrix with BlockLinearOperator in projectors and remove…

f9a92b8

… duplicated __del__ in StencilVector

spossann and others added 5 commits September 13, 2024 14:48

Update psydac/feec/tests/test_commuting_projections.py

a6ce393

Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

Fixes in test_3d_commuting_pro_1 (#438)

214079d

Change `IdentityLinearOperator` to `IdentityOperator`, and `e0 @ e0` to `e0.dot(e0)`. Add missing import statement.

Updated GlobalProjector so the m>1 is handled correctly

245a45e

Fixed bug, replaced imat with hmat

e66f863

yguclu assigned maxlin-ipp Sep 25, 2024

yguclu and others added 3 commits October 1, 2024 15:13

Improve and test Pyccel compiler flags (#428)

28e0f89

- 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

Update calls to numpy.reshape to avoid deprecation warnings (#433)

dd26dd5

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.

Circumvent apt-get download error (#440)

a3df57d

Update `apt` index files with `apt-get update` to avoid download error in CI workflow.

yguclu and others added 6 commits October 1, 2024 15:13

Remove is_block property from VectorFemSpace (#439)

c617b15

The obsolete property `is_block` was never used in Psydac nor in Struphy, hence we remove it. --------- Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

Add sqrt option to diagonal method of stencil/block matrices (#434)

2dc2297

--------- Co-authored-by: Yaman Güçlü <yaman.guclu@gmail.com>

Merge remote-tracking branch 'origin/devel' into expose-inter-histopo…

08b003d

…lation-matrices

reverted stencil2coo_kernels.py to devel version

6783d3f

yguclu requested changes Oct 2, 2024

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

spossann commented Oct 12, 2023 •

edited

Loading

spossann commented Oct 26, 2023

spossann commented Nov 20, 2023

yguclu commented Dec 5, 2023

spossann commented Dec 15, 2023

yguclu Oct 2, 2024

yguclu Oct 2, 2024

yguclu Oct 2, 2024

yguclu Oct 2, 2024

yguclu Oct 2, 2024

yguclu Oct 2, 2024

yguclu Oct 2, 2024

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

Are you sure you want to change the base?

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

Conversation

spossann commented Oct 12, 2023 • edited Loading

spossann commented Oct 26, 2023

spossann commented Nov 20, 2023

yguclu commented Dec 5, 2023

spossann commented Dec 15, 2023

yguclu Oct 2, 2024

Choose a reason for hiding this comment

yguclu Oct 2, 2024

Choose a reason for hiding this comment

yguclu Oct 2, 2024

Choose a reason for hiding this comment

yguclu Oct 2, 2024

Choose a reason for hiding this comment

yguclu Oct 2, 2024

Choose a reason for hiding this comment

yguclu Oct 2, 2024

Choose a reason for hiding this comment

yguclu Oct 2, 2024

Choose a reason for hiding this comment

spossann commented Oct 12, 2023 •

edited

Loading