From a6d4f998cd5b7ee4e98ec60c47205a53e9b2dc48 Mon Sep 17 00:00:00 2001 From: elmosa Date: Wed, 30 Apr 2025 20:50:00 +0200 Subject: [PATCH 01/25] fix issue of real scalar and complex domain --- psydac/linalg/basic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index 80bcac3d0..51101c458 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -546,7 +546,7 @@ def __init__(self, domain, codomain, c, A): assert isinstance(domain, VectorSpace) assert isinstance(codomain, VectorSpace) assert np.isscalar(c) - assert np.iscomplexobj(c) == (codomain._dtype == complex) + assert not np.iscomplexobj(c) or codomain._dtype == complex assert isinstance(A, LinearOperator) assert domain == A.domain assert codomain == A.codomain From 61760524e7a3e983be5da1373d272b5741d5f46c Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 5 May 2025 12:04:39 +0200 Subject: [PATCH 02/25] allow out argument in array_to_psydac --- psydac/linalg/utilities.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/utilities.py b/psydac/linalg/utilities.py index 57a9a6b86..eefad0e7b 100644 --- a/psydac/linalg/utilities.py +++ b/psydac/linalg/utilities.py @@ -15,7 +15,7 @@ ) #============================================================================== -def array_to_psydac(x, V): +def array_to_psydac(x, V, out=None): """ Convert a NumPy array to a Vector of the space V. This function is designed to be the inverse of the method .toarray() of the class Vector. Note: This function works in parallel but it is very costly and should be avoided if performance is a priority. @@ -28,6 +28,9 @@ def array_to_psydac(x, V): V : psydac.linalg.stencil.StencilVectorSpace or psydac.linalg.block.BlockVectorSpace Space of the final Psydac Vector. + out : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector, optional + The Psydac vector where to store the result. + Returns ------- u : psydac.linalg.stencil.StencilVector or psydac.linalg.block.BlockVector @@ -40,7 +43,13 @@ def array_to_psydac(x, V): assert V.dtype==complex, 'Complex array cannot be converted to a real StencilVector' assert x.size == V.dimension, 'Array must have the same global size as the space.' - u = V.zeros() + if out is None: + u = V.zeros() + else: + assert isinstance(out, BlockVector) + assert out.space is V + u = out + _array_to_psydac_recursive(x, u) u.update_ghost_regions() From 2b23fc09d1b6cb69dee27c91e8d27e349db7d93d Mon Sep 17 00:00:00 2001 From: elmosa Date: Sun, 11 May 2025 11:35:40 +0200 Subject: [PATCH 03/25] pass dtype to direct solver for interpolation and histopolation --- psydac/fem/splines.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psydac/fem/splines.py b/psydac/fem/splines.py index 03a1703cd..d6d02f593 100644 --- a/psydac/fem/splines.py +++ b/psydac/fem/splines.py @@ -182,7 +182,7 @@ def init_interpolation( self, dtype=float ): if self.periodic: # Convert to CSC format and compute sparse LU decomposition - self._interpolator = SparseSolver( csc_matrix( imat ) ) + self._interpolator = SparseSolver( csc_matrix(imat, dtype=dtype) ) else: # Convert to LAPACK banded format (see DGBTRF function) dmat = dia_matrix( imat ) @@ -218,7 +218,7 @@ def init_histopolation( self, dtype=float): self.hmat= imat if self.periodic: # Convert to CSC format and compute sparse LU decomposition - self._histopolator = SparseSolver( csc_matrix( imat ) ) + self._histopolator = SparseSolver( csc_matrix(imat, dtype=dtype) ) else: # Convert to LAPACK banded format (see DGBTRF function) dmat = dia_matrix( imat ) From bdd9c2e727e89ab5339dad0c5c9282ee461cce51 Mon Sep 17 00:00:00 2001 From: elmosa Date: Sat, 8 Nov 2025 18:55:00 +0100 Subject: [PATCH 04/25] improve to petsc --- psydac/linalg/block.py | 4 +- psydac/linalg/stencil.py | 4 +- psydac/linalg/tests/test_block.py | 31 +++++++++--- psydac/linalg/tests/test_stencil_vector.py | 58 ++++++++++++++++------ psydac/linalg/topetsc.py | 34 ++++++++----- 5 files changed, 93 insertions(+), 38 deletions(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index f2a755015..e2739b25d 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -481,11 +481,11 @@ def toarray_local(self, order='C'): return np.block([blocks])[0] # ... - def topetsc(self): + def topetsc(self, out=None): """ Convert to petsc data structure. """ from psydac.linalg.topetsc import vec_topetsc - vec = vec_topetsc( self ) + vec = vec_topetsc( self, out=out ) return vec #=============================================================================== diff --git a/psydac/linalg/stencil.py b/psydac/linalg/stencil.py index 73599655a..11e03aa55 100644 --- a/psydac/linalg/stencil.py +++ b/psydac/linalg/stencil.py @@ -698,11 +698,11 @@ def _toarray_parallel_with_pads(self, order='C'): return out.flatten( order=order) #... - def topetsc(self): + def topetsc(self, out=None): """ Convert to petsc data structure. """ from psydac.linalg.topetsc import vec_topetsc - vec = vec_topetsc( self ) + vec = vec_topetsc( self, out=out ) return vec # ... diff --git a/psydac/linalg/tests/test_block.py b/psydac/linalg/tests/test_block.py index ff3bf1e6c..327e6a40f 100644 --- a/psydac/linalg/tests/test_block.py +++ b/psydac/linalg/tests/test_block.py @@ -813,6 +813,8 @@ def test_block_vector_2d_serial_topetsc( dtype, n1, n2, p1, p2, P1, P2 ): W = BlockVectorSpace(V1, V2) x = BlockVector(W) + y = BlockVector(W) + z = BlockVector(W) # Fill in vector with random values, then update ghost regions for i1 in range(n1): @@ -823,10 +825,16 @@ def test_block_vector_2d_serial_topetsc( dtype, n1, n2, p1, p2, P1, P2 ): x.update_ghost_regions() v = x.topetsc() - v = petsc_to_psydac(v, W) + u = v.copy() + x.topetsc(out=u) - # The vectors can only be compared in the serial case - assert np.allclose( x.toarray() , v.toarray() ) + # Convert PETSc.Vec to BlockVector of W + y = petsc_to_psydac(v, W) + petsc_to_psydac(u, W, out=z) + + assert np.array_equal(x.toarray(), y.toarray()) + assert np.array_equal(x.toarray(), z.toarray()) + assert np.array_equal(v.array, u.array) #=============================================================================== @pytest.mark.parametrize( 'dtype', [float, complex] ) @@ -1245,6 +1253,8 @@ def test_block_vector_2d_parallel_topetsc( dtype, n1, n2, p1, p2, P1, P2 ): #TODO: implement conversion to PETSc recursively to treat case of block of blocks x = BlockVector(W) + y = BlockVector(W) + z = BlockVector(W) # Fill in vector with random values, then update ghost regions for i0 in range(len(W.starts)): @@ -1254,9 +1264,18 @@ def test_block_vector_2d_parallel_topetsc( dtype, n1, n2, p1, p2, P1, P2 ): x.update_ghost_regions() - v = petsc_to_psydac(x.topetsc(), W) + v = x.topetsc() + u = v.copy() + x.topetsc(out=u) + + # Convert PETSc.Vec to BlockVector of W + y = petsc_to_psydac(v, W) + petsc_to_psydac(u, W, out=z) + + assert np.allclose( x.toarray() , y.toarray(), rtol=1e-12, atol=1e-12 ) + assert np.allclose( x.toarray() , z.toarray(), rtol=1e-12, atol=1e-12 ) + assert np.allclose( v.array, u.array, rtol=1e-12, atol=1e-12 ) - assert np.allclose( x.toarray() , v.toarray(), rtol=1e-12, atol=1e-12 ) #=============================================================================== @pytest.mark.parametrize( 'dtype', [float, complex] ) @@ -1510,4 +1529,4 @@ def test_block_matrix_operator_parallel_dot_backend( dtype, n1, n2, p1, p2, P1, #=============================================================================== if __name__ == "__main__": import sys - pytest.main( sys.argv ) + pytest.main( sys.argv ) \ No newline at end of file diff --git a/psydac/linalg/tests/test_stencil_vector.py b/psydac/linalg/tests/test_stencil_vector.py index 0983edf80..d759ab106 100644 --- a/psydac/linalg/tests/test_stencil_vector.py +++ b/psydac/linalg/tests/test_stencil_vector.py @@ -429,6 +429,8 @@ def test_stencil_vector_2d_serial_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P2) # Create vector space and stencil vector V = StencilVectorSpace(C, dtype=dtype) x = StencilVector(V) + y = StencilVector(V) + z = StencilVector(V) # Fill the vector with data if dtype == complex: @@ -441,19 +443,24 @@ def test_stencil_vector_2d_serial_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P2) # Convert vector to PETSc.Vec v = x.topetsc() + w = v.copy() + x.topetsc(out=w) # Convert PETSc.Vec to StencilVector of V - v = petsc_to_psydac(v, V) + z = petsc_to_psydac(v, V) + petsc_to_psydac(w, V, out=y) # Test properties of v and data contained - assert v.space is V - assert v.dtype == dtype - assert v.starts == (0, 0) - assert v.ends == (n1 - 1, n2 - 1) - assert v.pads == (p1, p2) - assert v._data.shape == (n1 + 2 * p1 * s1, n2 + 2 * p2 * s2) - assert v._data.dtype == dtype - assert np.array_equal(x.toarray(), v.toarray()) + assert z.space is V + assert z.dtype == dtype + assert z.starts == (0, 0) + assert z.ends == (n1 - 1, n2 - 1) + assert z.pads == (p1, p2) + assert z._data.shape == (n1 + 2 * p1 * s1, n2 + 2 * p2 * s2) + assert z._data.dtype == dtype + assert np.array_equal(x.toarray(), z.toarray()) + assert np.array_equal(x.toarray(), y.toarray()) + assert np.array_equal(v.array, w.array) # =============================================================================== @pytest.mark.parametrize('dtype', [float, complex]) @@ -616,6 +623,8 @@ def test_stencil_vector_2d_parallel_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P # Create vector space and stencil vector V = StencilVectorSpace(C, dtype=dtype) x = StencilVector(V) + y = StencilVector(V) + z = StencilVector(V) # Fill the vector with data if dtype == complex: @@ -630,11 +639,16 @@ def test_stencil_vector_2d_parallel_topetsc(dtype, n1, n2, p1, p2, s1, s2, P1, P # Convert vector to PETSc.Vec v = x.topetsc() + u = v.copy() + x.topetsc(out=u) # Convert PETSc.Vec to StencilVector of V - v = petsc_to_psydac(v, V) + y = petsc_to_psydac(v, V) + petsc_to_psydac(u, V, out=z) - assert np.array_equal(x.toarray(), v.toarray()) + assert np.array_equal(x.toarray(), y.toarray()) + assert np.array_equal(x.toarray(), z.toarray()) + assert np.array_equal(v.array, u.array) # =============================================================================== @pytest.mark.parametrize('dtype', [float, complex]) @@ -659,6 +673,8 @@ def test_stencil_vector_1d_parallel_topetsc(dtype, n1, p1, s1, P1): # Create vector space and stencil vector V = StencilVectorSpace(C, dtype=dtype) x = StencilVector(V) + y = StencilVector(V) + z = StencilVector(V) # Fill the vector with data if dtype == complex: @@ -672,11 +688,16 @@ def test_stencil_vector_1d_parallel_topetsc(dtype, n1, p1, s1, P1): # Convert vector to PETSc.Vec v = x.topetsc() + u = v.copy() + x.topetsc(out=u) # Convert PETSc.Vec to StencilVector of V - v = petsc_to_psydac(v, V) + y = petsc_to_psydac(v, V) + petsc_to_psydac(u, V, out=z) - assert np.array_equal(x.toarray(), v.toarray()) + assert np.array_equal(x.toarray(), y.toarray()) + assert np.array_equal(x.toarray(), z.toarray()) + assert np.array_equal(v.array, u.array) # =============================================================================== @pytest.mark.parametrize('dtype', [float, complex]) @@ -710,6 +731,8 @@ def test_stencil_vector_3d_parallel_topetsc(dtype, n1, n2, n3, p1, p2, p3, s1, s # Create vector space and stencil vector V = StencilVectorSpace(C, dtype=dtype) x = StencilVector(V) + y = StencilVector(V) + z = StencilVector(V) # Fill the vector with data @@ -726,11 +749,16 @@ def test_stencil_vector_3d_parallel_topetsc(dtype, n1, n2, n3, p1, p2, p3, s1, s # Convert vector to PETSc.Vec v = x.topetsc() + u = v.copy() + x.topetsc(out=u) # Convert PETSc.Vec to StencilVector of V - v = petsc_to_psydac(v, V) + y = petsc_to_psydac(v, V) + petsc_to_psydac(u, V, out=z) - assert np.array_equal(x.toarray(), v.toarray()) + assert np.array_equal(x.toarray(), y.toarray()) + assert np.array_equal(x.toarray(), z.toarray()) + assert np.array_equal(v.array, u.array) # =============================================================================== @pytest.mark.parametrize('dtype', [float, complex]) diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py index d31cb5ca8..4b9269c1c 100644 --- a/psydac/linalg/topetsc.py +++ b/psydac/linalg/topetsc.py @@ -2,7 +2,7 @@ import numpy as np -from psydac.linalg.basic import VectorSpace +from psydac.linalg.basic import VectorSpace, ZeroOperator from psydac.linalg.block import BlockVectorSpace, BlockVector, BlockLinearOperator from psydac.linalg.stencil import StencilVectorSpace, StencilVector, StencilMatrix from psydac.linalg.kernels.stencil2IJV_kernels import stencil2IJV_1d_C, stencil2IJV_2d_C, stencil2IJV_3d_C @@ -332,7 +332,7 @@ def get_npts_per_block(V : VectorSpace) -> list: return npts_local_per_block -def vec_topetsc(vec): +def vec_topetsc(vec, out=None): """ Convert vector from Psydac format to a PETSc.Vec object. Parameters @@ -370,13 +370,19 @@ def vec_topetsc(vec): # Sum over the blocks to get the total local size localsize = np.sum(np.prod(npts_local, axis=1)) - gvec = PETSc.Vec().create(comm=carts[0].global_comm) + if out is not None: + assert isinstance(out, PETSc.Vec) + assert out.getSizes() == (localsize, globalsize) + out.zeroEntries() + gvec = out + else: + gvec = PETSc.Vec().create(comm=carts[0].global_comm) - # Set global and local size: - gvec.setSizes(size=(localsize, globalsize)) + # Set global and local size: + gvec.setSizes(size=(localsize, globalsize)) - gvec.setFromOptions() - gvec.setUp() + gvec.setFromOptions() + gvec.setUp() petsc_indices = [] petsc_data = [] @@ -506,12 +512,14 @@ def mat_topetsc(mat): for bc, bd in nonzero_block_indices: if isinstance(mat, BlockLinearOperator): mat_block = mat.blocks[bc][bd] - dnpts_block = dnpts_per_block_per_process[bd] - cnpts_block = cnpts_per_block_per_process[bc] - dshift_block = dindex_shift[bd] - cshift_block = cindex_shift[bc] + + if not isinstance(mat_block, ZeroOperator): + dnpts_block = dnpts_per_block_per_process[bd] + cnpts_block = cnpts_per_block_per_process[bc] + dshift_block = dindex_shift[bd] + cshift_block = cindex_shift[bc] - I,J,V,rowmap = toIJVrowmap(mat_block, bd, bc, I, J, V, rowmap, mat.domain, mat.codomain, dnpts_block, cnpts_block, dshift_block, cshift_block) + I,J,V,rowmap = toIJVrowmap(mat_block, bd, bc, I, J, V, rowmap, mat.domain, mat.codomain, dnpts_block, cnpts_block, dshift_block, cshift_block) # Set the values using IJV&rowmap format. The values are stored in a cache memory. gmat.setValuesIJV(I, J, V, rowmap=rowmap, addv=PETSc.InsertMode.ADD_VALUES) # The addition mode is necessary when periodic BC @@ -519,4 +527,4 @@ def mat_topetsc(mat): # Assemble the matrix with the values from the cache. Here it is where PETSc exchanges global communication. gmat.assemble() - return gmat + return gmat \ No newline at end of file From 2e69661e11534c737a7ac02a8b41f46ffaffe9d7 Mon Sep 17 00:00:00 2001 From: elmosa Date: Sat, 8 Nov 2025 20:13:27 +0100 Subject: [PATCH 05/25] avoid calling update_ghost_regions, exchange_assembly_data, remove_spurious_entries when the method does not exist --- psydac/linalg/block.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index e2739b25d..6dcc6a184 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -885,17 +885,20 @@ def nonzero_block_indices(self): # ... def update_ghost_regions(self): for Lij in self._blocks.values(): - Lij.update_ghost_regions() + if hasattr(Lij, 'update_ghost_regions'): + Lij.update_ghost_regions() # ... def exchange_assembly_data(self): for Lij in self._blocks.values(): - Lij.exchange_assembly_data() + if hasattr(Lij, 'exchange_assembly_data'): + Lij.exchange_assembly_data() # ... def remove_spurious_entries(self ): for Lij in self._blocks.values(): - Lij.remove_spurious_entries() + if hasattr(Lij, 'remove_spurious_entries'): + Lij.remove_spurious_entries() @property def ghost_regions_in_sync(self): From 7b96c83d2711e6c1172a170017b9bb2fa212d3e5 Mon Sep 17 00:00:00 2001 From: elmosa Date: Sat, 8 Nov 2025 20:14:51 +0100 Subject: [PATCH 06/25] delete space --- psydac/linalg/block.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/block.py b/psydac/linalg/block.py index 6dcc6a184..95db2d885 100644 --- a/psydac/linalg/block.py +++ b/psydac/linalg/block.py @@ -895,7 +895,7 @@ def exchange_assembly_data(self): Lij.exchange_assembly_data() # ... - def remove_spurious_entries(self ): + def remove_spurious_entries(self): for Lij in self._blocks.values(): if hasattr(Lij, 'remove_spurious_entries'): Lij.remove_spurious_entries() From 7c8345e96eeca883f2b6c21cc534e89dbdb8083d Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 7 Sep 2026 15:24:27 +0200 Subject: [PATCH 07/25] Avoid running into max_it when a location coincides exactly with the right bound of the domain --- psydac/core/bsplines_kernels.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index caaf75a95..5d6d94962 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -1095,7 +1095,7 @@ def cell_index_p(breaks: 'float[:]', i_grid: 'float[:]', tol: float, out: 'int[: i_cell = (low + high)//2 it = 0 max_it = 2*(high-low) # this number of iterations should not be reached - while it < max_it and (x < breaks[i_cell] - tol or x >= breaks[i_cell + 1] + tol): + while it < max_it and (x < breaks[i_cell] - tol or x > breaks[i_cell + 1] + tol): if x < breaks[i_cell]: high = i_cell else: @@ -1103,6 +1103,7 @@ def cell_index_p(breaks: 'float[:]', i_grid: 'float[:]', tol: float, out: 'int[: i_cell = (low + high)//2 it += 1 if it >= max_it: + print(current_index, x, i_cell) return -2 # Check were we landed with the binary search From 53423194f7500070d44127a917d4f05db0838fa4 Mon Sep 17 00:00:00 2001 From: elmosa Date: Tue, 8 Sep 2026 16:41:22 +0200 Subject: [PATCH 08/25] delete forgotten print --- psydac/core/bsplines_kernels.py | 1 - 1 file changed, 1 deletion(-) diff --git a/psydac/core/bsplines_kernels.py b/psydac/core/bsplines_kernels.py index 5d6d94962..105bf1a95 100644 --- a/psydac/core/bsplines_kernels.py +++ b/psydac/core/bsplines_kernels.py @@ -1103,7 +1103,6 @@ def cell_index_p(breaks: 'float[:]', i_grid: 'float[:]', tol: float, out: 'int[: i_cell = (low + high)//2 it += 1 if it >= max_it: - print(current_index, x, i_cell) return -2 # Check were we landed with the binary search From 666024cc530c82d4a5eb82b322d9b9e98439013b Mon Sep 17 00:00:00 2001 From: elmosa Date: Fri, 18 Sep 2026 15:26:32 +0200 Subject: [PATCH 09/25] add topetsc method for ZeroOperator, and a test --- psydac/linalg/basic.py | 6 ++++- psydac/linalg/tests/test_matrix_free.py | 25 ++++++++++++++++++++ psydac/linalg/topetsc.py | 31 ++++++++++++++----------- 3 files changed, 48 insertions(+), 14 deletions(-) diff --git a/psydac/linalg/basic.py b/psydac/linalg/basic.py index a6bc99b66..402afc67f 100644 --- a/psydac/linalg/basic.py +++ b/psydac/linalg/basic.py @@ -553,12 +553,16 @@ def copy(self): return ZeroOperator(self.domain, self.codomain) def toarray(self): - return np.zeros(self.shape, dtype=self.dtype) + return np.zeros(self.shape, dtype=self.dtype) def tosparse(self): from scipy.sparse import csr_matrix return csr_matrix(self.shape, dtype=self.dtype) + def topetsc(self): + from psydac.linalg.topetsc import mat_topetsc + return mat_topetsc(self) + def transpose(self, conjugate=False): return ZeroOperator(domain=self.codomain, codomain=self.domain) diff --git a/psydac/linalg/tests/test_matrix_free.py b/psydac/linalg/tests/test_matrix_free.py index 1bf1a9fa7..3e28e36d9 100644 --- a/psydac/linalg/tests/test_matrix_free.py +++ b/psydac/linalg/tests/test_matrix_free.py @@ -128,6 +128,31 @@ def test_solvers_matrix_free(solver, use_jacobi_pc): error = np.linalg.norm( (b - y).toarray()) assert np.linalg.norm( (b - y).toarray() ) < tol +@pytest.mark.parametrize('n1', [3, 5]) +@pytest.mark.parametrize('n2', [4, 7]) +@pytest.mark.parametrize('p1', [2, 6]) +@pytest.mark.parametrize('p2', [3, 9]) +@pytest.mark.petsc +def test_ZeroOperator_topetsc(n1, n2, p1, p2): + P1 = False + P2 = False + m1 = (n2+n1)//2 + m2 = n1+1 + q1 = p1 + q2 = p2 + V1 = get_StencilVectorSpace([n1, n2], [p1, p2], [P1, P2]) + V2 = get_StencilVectorSpace([m1, m2], [q1, q2], [P1, P2]) + O = ZeroOperator(codomain=V2, domain=V1) + Op = O.topetsc() + + v1 = get_random_StencilVector(V1) + v1p = v1.topetsc() + v2p = Op.createVecLeft() + + Op.mult(v1p, v2p) + + assert np.linalg.norm( v2p.array ) == 0. + #=============================================================================== # SCRIPT FUNCTIONALITY #=============================================================================== diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py index 666b8a682..8865c963a 100644 --- a/psydac/linalg/topetsc.py +++ b/psydac/linalg/topetsc.py @@ -458,7 +458,7 @@ def mat_topetsc(mat): from petsc4py import PETSc - assert isinstance(mat, StencilMatrix) or isinstance(mat, BlockLinearOperator), 'Conversion only implemented for StencilMatrix and BlockLinearOperator.' + assert isinstance(mat, StencilMatrix) or isinstance(mat, BlockLinearOperator) or isinstance(mat, ZeroOperator), 'Conversion only implemented for StencilMatrix, BlockLinearOperator and ZeroOperator.' if (isinstance(mat.domain, BlockVectorSpace) and any([isinstance(mat.domain.spaces[b], BlockVectorSpace) for b in range(len(mat.domain.spaces))]))\ or (isinstance(mat.codomain, BlockVectorSpace) and any([isinstance(mat.codomain.spaces[b], BlockVectorSpace) for b in range(len(mat.codomain.spaces))])): @@ -469,23 +469,11 @@ def mat_topetsc(mat): elif isinstance(mat.domain, BlockVectorSpace): comm = mat.domain.spaces[0].cart.global_comm - nonzero_block_indices = ((0,0),) if isinstance(mat, StencilMatrix) else mat.nonzero_block_indices - - mat.update_ghost_regions() - mat.remove_spurious_entries() # Get the number of points local to the current process: dnpts_local = get_npts_local(mat.domain) # indexed [block, dimension]. Different for each process. cnpts_local = get_npts_local(mat.codomain) # indexed [block, dimension]. Different for each process. - # Get the number of points per block, per process and per dimension: - dnpts_per_block_per_process = np.array(get_npts_per_block(mat.domain)) # global variable, indexed as [block, process, dimension] - cnpts_per_block_per_process = np.array(get_npts_per_block(mat.codomain)) # global variable, indexed as [block, process, dimension] - - # Get the index shift for each block and each process: - dindex_shift = get_index_shift_per_block_per_process(mat.domain) # global variable, indexed as [block, process, dimension] - cindex_shift = get_index_shift_per_block_per_process(mat.codomain) # global variable, indexed as [block, process, dimension] - globalsize = mat.shape # Sum over the blocks to get the total local size @@ -506,6 +494,23 @@ def mat_topetsc(mat): gmat.setFromOptions() gmat.setUp() + if isinstance(mat, ZeroOperator): + gmat.assemble() + return gmat + + nonzero_block_indices = ((0,0),) if isinstance(mat, StencilMatrix) else mat.nonzero_block_indices + + mat.update_ghost_regions() + mat.remove_spurious_entries() + + # Get the number of points per block, per process and per dimension: + dnpts_per_block_per_process = np.array(get_npts_per_block(mat.domain)) # global variable, indexed as [block, process, dimension] + cnpts_per_block_per_process = np.array(get_npts_per_block(mat.codomain)) # global variable, indexed as [block, process, dimension] + + # Get the index shift for each block and each process: + dindex_shift = get_index_shift_per_block_per_process(mat.domain) # global variable, indexed as [block, process, dimension] + cindex_shift = get_index_shift_per_block_per_process(mat.codomain) # global variable, indexed as [block, process, dimension] + I = [0] # Row pointers J = [] # Column indices V = [] # Values From dfb0fb6512a4d8341adaa3fd4b4c1f51712ac188 Mon Sep 17 00:00:00 2001 From: elmosa Date: Fri, 18 Sep 2026 16:46:12 +0200 Subject: [PATCH 10/25] test projectors with complex dtype --- psydac/feec/tests/test_global_projectors.py | 32 +++++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index f4a549e38..6cd353654 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -39,13 +39,13 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') - V0 = TensorFemSpace(domain_decomposition, N) + V0 = TensorFemSpace(domain_decomposition, N, dtype=complex) # Projector onto H1 space (1D interpolation) P0 = GlobalGeometricProjectorH1(V0) # Function to project - f = lambda xi1 : np.sin( xi1 + 0.5 ) + f = lambda xi1 : np.sin( (2+1j)*xi1 + 0.5 ) # Compute the projection u0 = P0(f) @@ -86,7 +86,7 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, multiplicity = min(multiplicity, degree) # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') - V0 = TensorFemSpace(domain_decomposition, N) + V0 = TensorFemSpace(domain_decomposition, N, dtype=complex) # L2 space (1-forms) V1 = V0.reduce_degree(axes=[0], basis='M') @@ -95,7 +95,7 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, P1 = GlobalGeometricProjectorL2(V1, nquads=[nquads]) # Function to project - f = lambda xi1 : np.sin( xi1 + 0.5 ) + f = lambda xi1 : np.sin( xi1 + 0.5j ) # Compute the projection u1 = P1(f) @@ -129,11 +129,15 @@ def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbos multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)] derham = Derham(domain, ["H1", "Hdiv", "L2"]) + derham.V0.codomain_type='complex' + derham.V1.codomain_type='complex' + derham.V2.codomain_type='complex' + derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity) P0, P1, P2, PX = derham_h.projectors(nquads=[2*p+1 for p in degree]) # Function to project - f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) + f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( (1+1e-2j)*xi2 + 0.3 ) f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) # Compute the projection @@ -193,11 +197,14 @@ def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, verb multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)] derham = Derham(domain, ["H1", "Hdiv", "L2"]) + derham.V0.codomain_type='complex' + derham.V1.codomain_type='complex' + derham.V2.codomain_type='complex' derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity) P0, P1, P2, PX = derham_h.projectors() # Function to project - f1 = lambda xi1, xi2 : 20 * xi1**2*(xi1-1.)**2 + f1 = lambda xi1, xi2 : (20+2j) * xi1**2*(xi1-1.)**2 f2 = lambda xi1, xi2 : 10 * xi2**2*(xi2-1.)**2 # Compute the projection @@ -265,12 +272,15 @@ def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, verbo multiplicity = [min(m, p) for p, m in zip (degree, multiplicity)] derham = Derham(domain, ["H1", "Hcurl", "L2"]) + derham.V0.codomain_type='complex' + derham.V1.codomain_type='complex' + derham.V2.codomain_type='complex' derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity) P0, P1, P2, PX = derham_h.projectors() # Function to project - f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) - f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) + f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3j ) + f2 = lambda xi1, xi2 : 1j * np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) # Compute the projection u0 = P0(f1) @@ -324,6 +334,10 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal domain_h = discretize(domain, ncells=ncells, periodic=periodic) derham = Derham(domain) + derham.V0.codomain_type='complex' + derham.V1.codomain_type='complex' + derham.V2.codomain_type='complex' + derham.V3.codomain_type='complex' #change multiplicity if higher than degree to avoid problems (case p Date: Fri, 18 Sep 2026 19:07:56 +0200 Subject: [PATCH 11/25] do not use complex dtype in all the tests --- psydac/feec/tests/test_global_projectors.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index 6cd353654..07cb1f81e 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -39,13 +39,13 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') - V0 = TensorFemSpace(domain_decomposition, N, dtype=complex) + V0 = TensorFemSpace(domain_decomposition, N) # Projector onto H1 space (1D interpolation) P0 = GlobalGeometricProjectorH1(V0) # Function to project - f = lambda xi1 : np.sin( (2+1j)*xi1 + 0.5 ) + f = lambda xi1 : np.sin( xi1 + 0.5 ) # Compute the projection u0 = P0(f) @@ -129,15 +129,12 @@ def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbos multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)] derham = Derham(domain, ["H1", "Hdiv", "L2"]) - derham.V0.codomain_type='complex' - derham.V1.codomain_type='complex' - derham.V2.codomain_type='complex' derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity) P0, P1, P2, PX = derham_h.projectors(nquads=[2*p+1 for p in degree]) # Function to project - f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( (1+1e-2j)*xi2 + 0.3 ) + f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) # Compute the projection @@ -347,7 +344,7 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal # Function to project f1 = lambda xi1, xi2, xi3 : np.sin( xi1 + 0.51 ) * np.cos( xi2 + 0.32 ) * np.sin( xi3 - 0.43) f2 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.27 ) * np.sin( xi2 - 0.29 ) * np.cos( xi3 + 0.67) - f3 = lambda xi1, xi2, xi3 : np.cos( (1+1e-2j)*xi1 + 0.72 ) * np.sin( xi2 - 0.73j ) * np.cos( xi3 - 0.14) + f3 = lambda xi1, xi2, xi3 : np.cos( (1+1e-2j)*xi1 + 0.72 ) * np.sin( xi2 - 0.73 + 0.1j ) * np.cos( xi3 - 0.14) # Compute the projection u0 = P0(f1) From ef36b01dcefa9a32b54bea66c17ea263913d7275 Mon Sep 17 00:00:00 2001 From: elmosa Date: Fri, 18 Sep 2026 19:13:35 +0200 Subject: [PATCH 12/25] use periodic fun for periodic tests --- psydac/feec/tests/test_global_projectors.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index 07cb1f81e..643fc6c99 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -39,13 +39,13 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') - V0 = TensorFemSpace(domain_decomposition, N) + V0 = TensorFemSpace(domain_decomposition, N, dtype=complex) # Projector onto H1 space (1D interpolation) P0 = GlobalGeometricProjectorH1(V0) # Function to project - f = lambda xi1 : np.sin( xi1 + 0.5 ) + f = lambda xi1 : np.sin( xi1 + 0.5 -1j ) # Compute the projection u0 = P0(f) @@ -61,6 +61,7 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose if verbose: print(ncells, maxnorm_error / error_estim) assert maxnorm_error <= max(10 * error_estim, 1e-13) +test_H1_projector_1d((0, 2*np.pi), 500, 1, True,1) #============================================================================== @pytest.mark.parametrize('domain', [(0, 2*np.pi)]) @@ -135,7 +136,7 @@ def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbos # Function to project f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) - f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) + f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2j ) # Compute the projection u0 = P0(f1) @@ -344,7 +345,7 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal # Function to project f1 = lambda xi1, xi2, xi3 : np.sin( xi1 + 0.51 ) * np.cos( xi2 + 0.32 ) * np.sin( xi3 - 0.43) f2 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.27 ) * np.sin( xi2 - 0.29 ) * np.cos( xi3 + 0.67) - f3 = lambda xi1, xi2, xi3 : np.cos( (1+1e-2j)*xi1 + 0.72 ) * np.sin( xi2 - 0.73 + 0.1j ) * np.cos( xi3 - 0.14) + f3 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.72 ) * np.sin( xi2 - 0.73 + 7j ) * np.cos( xi3 - 0.14) # Compute the projection u0 = P0(f1) From 1be62271915730c11cf61bd7d09c2a8437705bdf Mon Sep 17 00:00:00 2001 From: elmosa Date: Tue, 22 Sep 2026 10:37:21 +0200 Subject: [PATCH 13/25] fix tests --- psydac/feec/tests/test_global_projectors.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index 643fc6c99..a36299f2e 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -61,7 +61,6 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose if verbose: print(ncells, maxnorm_error / error_estim) assert maxnorm_error <= max(10 * error_estim, 1e-13) -test_H1_projector_1d((0, 2*np.pi), 500, 1, True,1) #============================================================================== @pytest.mark.parametrize('domain', [(0, 2*np.pi)]) @@ -345,7 +344,7 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal # Function to project f1 = lambda xi1, xi2, xi3 : np.sin( xi1 + 0.51 ) * np.cos( xi2 + 0.32 ) * np.sin( xi3 - 0.43) f2 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.27 ) * np.sin( xi2 - 0.29 ) * np.cos( xi3 + 0.67) - f3 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.72 ) * np.sin( xi2 - 0.73 + 7j ) * np.cos( xi3 - 0.14) + f3 = lambda xi1, xi2, xi3 : 1j * np.cos( xi1 + 0.72 ) * np.sin( xi2 - 0.73 ) * np.cos( xi3 - 0.14) # Compute the projection u0 = P0(f1) @@ -407,7 +406,7 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal if verbose: print(ncells, maxnorm_error / error_estim) assert maxnorm_error <= max(15 * error_estim, 1e-13) - +test_derham_projector_3d([10,9,12], [2,2,2], [True, True, True], [2,2,2], verbose=True) #============================================================================== def manual_convergence_tests(dim): """ From 6f40db4a6e35cd2f3bbe0db9dd6104a28e368b90 Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 10:39:46 +0200 Subject: [PATCH 14/25] new line --- psydac/linalg/tests/test_block.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/tests/test_block.py b/psydac/linalg/tests/test_block.py index 442518a4e..6390cfd83 100644 --- a/psydac/linalg/tests/test_block.py +++ b/psydac/linalg/tests/test_block.py @@ -1637,4 +1637,4 @@ def test_block_matrix_operator_parallel_dot_backend( dtype, n1, n2, p1, p2, P1, #=============================================================================== if __name__ == "__main__": import sys - pytest.main( sys.argv ) \ No newline at end of file + pytest.main( sys.argv ) From 8e84c6a707292fbaa4ae05c1f45a1771976e3646 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Elena=20Moral=20S=C3=A1nchez?= <88042165+e-moral-sanchez@users.noreply.github.com> Date: Mon, 12 Oct 2026 10:42:07 +0200 Subject: [PATCH 15/25] Update psydac/linalg/tests/test_matrix_free.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- psydac/linalg/tests/test_matrix_free.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/tests/test_matrix_free.py b/psydac/linalg/tests/test_matrix_free.py index 3e28e36d9..49a657904 100644 --- a/psydac/linalg/tests/test_matrix_free.py +++ b/psydac/linalg/tests/test_matrix_free.py @@ -142,7 +142,7 @@ def test_ZeroOperator_topetsc(n1, n2, p1, p2): q2 = p2 V1 = get_StencilVectorSpace([n1, n2], [p1, p2], [P1, P2]) V2 = get_StencilVectorSpace([m1, m2], [q1, q2], [P1, P2]) - O = ZeroOperator(codomain=V2, domain=V1) + O = ZeroOperator(domain=V1, codomain=V2) Op = O.topetsc() v1 = get_random_StencilVector(V1) From 1399a780c288ef910ab5b97dd657bc19d7915dd9 Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 10:44:48 +0200 Subject: [PATCH 16/25] add line --- psydac/linalg/topetsc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py index 8865c963a..a8e7be314 100644 --- a/psydac/linalg/topetsc.py +++ b/psydac/linalg/topetsc.py @@ -536,4 +536,4 @@ def mat_topetsc(mat): # Assemble the matrix with the values from the cache. Here it is where PETSc exchanges global communication. gmat.assemble() - return gmat \ No newline at end of file + return gmat From 3b24598df6598af44f13d3404fdd9f6e54545429 Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 11:25:07 +0200 Subject: [PATCH 17/25] document out --- psydac/linalg/topetsc.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py index a8e7be314..a0c4b752b 100644 --- a/psydac/linalg/topetsc.py +++ b/psydac/linalg/topetsc.py @@ -343,7 +343,9 @@ def vec_topetsc(vec, out=None): ---------- vec : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector PSYDAC StencilVector or BlockVector. In the case of a BlockVector, only the case where the blocks are StencilVector is implemented. - + out : PETSc.Vec + PETSc vector where to store the returned value. In case of None, a new vector is created. + Returns ------- gvec : PETSc.Vec From a832ec5e5d4c8c0214ff9a704c553581388b594c Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 11:39:43 +0200 Subject: [PATCH 18/25] parametrize dtype in 1d 2d tests --- psydac/feec/tests/test_global_projectors.py | 62 +++++++++++++++------ 1 file changed, 44 insertions(+), 18 deletions(-) diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index a36299f2e..6cf4bfcdc 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -25,8 +25,9 @@ @pytest.mark.parametrize('degree', [1, 2, 3, 4, 5, 6, 7]) @pytest.mark.parametrize('periodic', [False, True]) @pytest.mark.parametrize('multiplicity', [1, 2]) +@pytest.mark.parametrize('dtype', [float, complex]) -def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose=False): +def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, dtype, verbose=False): if verbose: print('\n1d_h1') @@ -39,13 +40,16 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') - V0 = TensorFemSpace(domain_decomposition, N, dtype=complex) + V0 = TensorFemSpace(domain_decomposition, N, dtype=dtype) # Projector onto H1 space (1D interpolation) P0 = GlobalGeometricProjectorH1(V0) # Function to project - f = lambda xi1 : np.sin( xi1 + 0.5 -1j ) + if dtype == complex: + f = lambda xi1 : np.sin( xi1 + 0.5 -1j ) + else: + f = lambda xi1 : np.sin( xi1 + 0.5) # Compute the projection u0 = P0(f) @@ -69,8 +73,9 @@ def test_H1_projector_1d(domain, ncells, degree, periodic, multiplicity, verbose @pytest.mark.parametrize('periodic', [True, False]) @pytest.mark.parametrize('nquads', [100, 120, 140, 160]) @pytest.mark.parametrize('multiplicity', [1, 2]) +@pytest.mark.parametrize('dtype', [float, complex]) -def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, verbose=False): +def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, dtype, verbose=False): if verbose: print('\n1d_l2') @@ -86,7 +91,7 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, multiplicity = min(multiplicity, degree) # H1 space (0-forms) N = SplineSpace(degree=degree, knots=knots, periodic=periodic, basis='B') - V0 = TensorFemSpace(domain_decomposition, N, dtype=complex) + V0 = TensorFemSpace(domain_decomposition, N, dtype=dtype) # L2 space (1-forms) V1 = V0.reduce_degree(axes=[0], basis='M') @@ -95,7 +100,10 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, P1 = GlobalGeometricProjectorL2(V1, nquads=[nquads]) # Function to project - f = lambda xi1 : np.sin( xi1 + 0.5j ) + if dtype == complex: + f = lambda xi1 : np.sin( xi1 + 0.5j ) + else: + f = lambda xi1 : np.sin( xi1 + 0.5 ) # Compute the projection u1 = P1(f) @@ -117,8 +125,9 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads, multiplicity, @pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]]) @pytest.mark.parametrize('periodic', [[False, False], [True, True]]) @pytest.mark.parametrize('multiplicity', [(1, 1), (2, 2)]) +@pytest.mark.parametrize('dtype', [float, complex]) -def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbose=False): +def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, dtype, verbose=False): if verbose: print('\n2d_hdiv') @@ -129,13 +138,20 @@ def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbos multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)] derham = Derham(domain, ["H1", "Hdiv", "L2"]) + if dtype == complex: + derham.V0.codomain_type='complex' + derham.V1.codomain_type='complex' + derham.V2.codomain_type='complex' derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity) P0, P1, P2, PX = derham_h.projectors(nquads=[2*p+1 for p in degree]) # Function to project f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) - f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2j ) + if dtype == complex: + f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2j ) + else: + f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) # Compute the projection u0 = P0(f1) @@ -182,8 +198,9 @@ def test_derham_projector_2d_hdiv(ncells, degree, periodic, multiplicity, verbos @pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]]) @pytest.mark.parametrize('periodic', [[False, False], [True, True]]) @pytest.mark.parametrize('multiplicity', [(1, 1), (2, 2)]) +@pytest.mark.parametrize('dtype', [float, complex]) -def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, verbose=False): +def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, dtype, verbose=False): if verbose: print('\n2d_hdiv_2') @@ -194,14 +211,18 @@ def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, verb multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)] derham = Derham(domain, ["H1", "Hdiv", "L2"]) - derham.V0.codomain_type='complex' - derham.V1.codomain_type='complex' - derham.V2.codomain_type='complex' + if dtype == complex: + derham.V0.codomain_type='complex' + derham.V1.codomain_type='complex' + derham.V2.codomain_type='complex' derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity) P0, P1, P2, PX = derham_h.projectors() # Function to project - f1 = lambda xi1, xi2 : (20+2j) * xi1**2*(xi1-1.)**2 + if dtype == complex: + f1 = lambda xi1, xi2 : (20+2j) * xi1**2*(xi1-1.)**2 + else: + f1 = lambda xi1, xi2 : 20 * xi1**2*(xi1-1.)**2 f2 = lambda xi1, xi2 : 10 * xi2**2*(xi2-1.)**2 # Compute the projection @@ -257,8 +278,9 @@ def test_derham_projector_2d_hdiv_2(ncells, degree, periodic, multiplicity, verb @pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]]) @pytest.mark.parametrize('periodic', [[False, False], [True, False] ,[True, True]]) @pytest.mark.parametrize('multiplicity', [[1,1],[2,2]]) +@pytest.mark.parametrize('dtype', [float, complex]) -def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, verbose=False): +def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, dtype, verbose=False): if verbose: print('\n2d_hcurl') @@ -269,14 +291,18 @@ def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, verbo multiplicity = [min(m, p) for p, m in zip (degree, multiplicity)] derham = Derham(domain, ["H1", "Hcurl", "L2"]) - derham.V0.codomain_type='complex' - derham.V1.codomain_type='complex' - derham.V2.codomain_type='complex' + if dtype == complex: + derham.V0.codomain_type='complex' + derham.V1.codomain_type='complex' + derham.V2.codomain_type='complex' derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True, multiplicity=multiplicity) P0, P1, P2, PX = derham_h.projectors() # Function to project - f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3j ) + if dtype == complex: + f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3j ) + else: + f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) f2 = lambda xi1, xi2 : 1j * np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) # Compute the projection From 77830c17c6c2fb9023978b5cc0c4ec042d79fbf2 Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 12:57:59 +0200 Subject: [PATCH 19/25] fix comments --- psydac/linalg/topetsc.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py index a0c4b752b..9ac676e32 100644 --- a/psydac/linalg/topetsc.py +++ b/psydac/linalg/topetsc.py @@ -505,13 +505,13 @@ def mat_topetsc(mat): mat.update_ghost_regions() mat.remove_spurious_entries() - # Get the number of points per block, per process and per dimension: - dnpts_per_block_per_process = np.array(get_npts_per_block(mat.domain)) # global variable, indexed as [block, process, dimension] - cnpts_per_block_per_process = np.array(get_npts_per_block(mat.codomain)) # global variable, indexed as [block, process, dimension] + # Get the number of points per block, per process and per dimension (global variable) indexed as [block, process, dimension] + dnpts_per_block_per_process = np.array(get_npts_per_block(mat.domain)) + cnpts_per_block_per_process = np.array(get_npts_per_block(mat.codomain)) - # Get the index shift for each block and each process: - dindex_shift = get_index_shift_per_block_per_process(mat.domain) # global variable, indexed as [block, process, dimension] - cindex_shift = get_index_shift_per_block_per_process(mat.codomain) # global variable, indexed as [block, process, dimension] + # Get the index shift per block, per process and per dimension (global variable) indexed as [block, process, dimension] + dindex_shift = get_index_shift_per_block_per_process(mat.domain) + cindex_shift = get_index_shift_per_block_per_process(mat.codomain) I = [0] # Row pointers J = [] # Column indices From a361a4b744e7079b3aca83fba78384e2ddabe931 Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 13:01:37 +0200 Subject: [PATCH 20/25] use array_equal test petsc --- psydac/linalg/tests/test_block.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/psydac/linalg/tests/test_block.py b/psydac/linalg/tests/test_block.py index 6390cfd83..9760396b6 100644 --- a/psydac/linalg/tests/test_block.py +++ b/psydac/linalg/tests/test_block.py @@ -1373,9 +1373,9 @@ def test_block_vector_2d_parallel_topetsc( dtype, n1, n2, p1, p2, P1, P2 ): y = petsc_to_psydac(v, W) petsc_to_psydac(u, W, out=z) - assert np.allclose( x.toarray() , y.toarray(), rtol=1e-12, atol=1e-12 ) - assert np.allclose( x.toarray() , z.toarray(), rtol=1e-12, atol=1e-12 ) - assert np.allclose( v.array, u.array, rtol=1e-12, atol=1e-12 ) + assert np.array_equal(x.toarray(), y.toarray()) + assert np.array_equal(x.toarray(), z.toarray()) + assert np.array_equal(v.array, u.array) #=============================================================================== From 4f0515618d02993f37c53b085c975e554b76b406 Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 13:05:14 +0200 Subject: [PATCH 21/25] skip one iteration if block is zero operator --- psydac/linalg/topetsc.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py index 9ac676e32..acb4d5c91 100644 --- a/psydac/linalg/topetsc.py +++ b/psydac/linalg/topetsc.py @@ -523,14 +523,16 @@ def mat_topetsc(mat): for bc, bd in nonzero_block_indices: if isinstance(mat, BlockLinearOperator): mat_block = mat.blocks[bc][bd] - - if not isinstance(mat_block, ZeroOperator): - dnpts_block = dnpts_per_block_per_process[bd] - cnpts_block = cnpts_per_block_per_process[bc] - dshift_block = dindex_shift[bd] - cshift_block = cindex_shift[bc] - I,J,V,rowmap = toIJVrowmap(mat_block, bd, bc, I, J, V, rowmap, mat.domain, mat.codomain, dnpts_block, cnpts_block, dshift_block, cshift_block) + if isinstance(mat_block, ZeroOperator): + continue + + dnpts_block = dnpts_per_block_per_process[bd] + cnpts_block = cnpts_per_block_per_process[bc] + dshift_block = dindex_shift[bd] + cshift_block = cindex_shift[bc] + + I,J,V,rowmap = toIJVrowmap(mat_block, bd, bc, I, J, V, rowmap, mat.domain, mat.codomain, dnpts_block, cnpts_block, dshift_block, cshift_block) # Set the values using IJV&rowmap format. The values are stored in a cache memory. gmat.setValuesIJV(I, J, V, rowmap=rowmap, addv=PETSc.InsertMode.ADD_VALUES) # The addition mode is necessary when periodic BC From a0a104e9c1aafbf98f6e33a4300914d089b88ce2 Mon Sep 17 00:00:00 2001 From: elmosa Date: Mon, 12 Oct 2026 17:32:08 +0200 Subject: [PATCH 22/25] forgot to put complex function within complex case --- psydac/feec/tests/test_global_projectors.py | 23 ++++++++++++--------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index 6cf4bfcdc..246179f36 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -301,9 +301,10 @@ def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, dtype # Function to project if dtype == complex: f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3j ) + f2 = lambda xi1, xi2 : 1j * np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) else: f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) - f2 = lambda xi1, xi2 : 1j * np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) + f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) # Compute the projection u0 = P0(f1) @@ -344,7 +345,7 @@ def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, dtype if verbose: print(ncells, maxnorm_error / error_estim) assert maxnorm_error <= max(10 * error_estim, 1e-13) - + #============================================================================== @pytest.mark.parametrize('ncells', [[10,9,12]]) @pytest.mark.parametrize('degree', [[2,2,2], [2,3,2], [3,3,3]]) @@ -432,7 +433,7 @@ def test_derham_projector_3d(ncells, degree, periodic, multiplicity, verbose=Fal if verbose: print(ncells, maxnorm_error / error_estim) assert maxnorm_error <= max(15 * error_estim, 1e-13) -test_derham_projector_3d([10,9,12], [2,2,2], [True, True, True], [2,2,2], verbose=True) + #============================================================================== def manual_convergence_tests(dim): """ @@ -444,23 +445,25 @@ def manual_convergence_tests(dim): if dim == 1: degrees = [3,4,5,6,7] ncells = [10, 20, 40, 80, 160, 320, 640] + dtype = float for deg in degrees: print(f"degree = {deg}") for nc in ncells: - test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=1, verbose=True) - test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=2, verbose=True) - test_L2_projector_1d(domain, nc, deg, periodic, nquads=deg, multiplicity=1, verbose=True) + test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=1, dtype=dtype, verbose=True) + test_H1_projector_1d(domain, nc, deg, periodic, multiplicity=2, dtype=dtype, verbose=True) + test_L2_projector_1d(domain, nc, deg, periodic, nquads=deg, multiplicity=1, dtype=dtype, verbose=True) elif dim == 2: degrees = [[3,2], [3, 3], [3,4]] ncells = [10, 20, 40, 80, 160] + dtype = float for deg in degrees: print(f"degree = {deg}") for nc in ncells: - test_derham_projector_2d_hcurl ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], verbose=True) - test_derham_projector_2d_hdiv ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], verbose=True) - test_derham_projector_2d_hdiv ([nc, nc], deg, [periodic, periodic], multiplicity=[2,2], verbose=True) - test_derham_projector_2d_hdiv_2([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], verbose=True) + test_derham_projector_2d_hcurl ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], dtype=dtype, verbose=True) + test_derham_projector_2d_hdiv ([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], dtype=dtype, verbose=True) + test_derham_projector_2d_hdiv ([nc, nc], deg, [periodic, periodic], multiplicity=[2,2], dtype=dtype, verbose=True) + test_derham_projector_2d_hdiv_2([nc, nc], deg, [periodic, periodic], multiplicity=[1,1], dtype=dtype, verbose=True) elif dim == 3: print("\n3d") From 4942adf421598c100d0381670764598fdfa76958 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Elena=20Moral=20S=C3=A1nchez?= <88042165+e-moral-sanchez@users.noreply.github.com> Date: Tue, 13 Oct 2026 11:35:06 +0200 Subject: [PATCH 23/25] Update psydac/linalg/topetsc.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Yaman Güçlü --- psydac/linalg/topetsc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/psydac/linalg/topetsc.py b/psydac/linalg/topetsc.py index acb4d5c91..44a9afa59 100644 --- a/psydac/linalg/topetsc.py +++ b/psydac/linalg/topetsc.py @@ -343,8 +343,8 @@ def vec_topetsc(vec, out=None): ---------- vec : psydac.linalg.stencil.StencilVector | psydac.linalg.block.BlockVector PSYDAC StencilVector or BlockVector. In the case of a BlockVector, only the case where the blocks are StencilVector is implemented. - out : PETSc.Vec - PETSc vector where to store the returned value. In case of None, a new vector is created. + out : PETSc.Vec, optional + PETSc vector where to store the returned value. In case of `None`, a new vector is created. Returns ------- From 08cf604b5f04e0e4499e9a8e6c82c4eb0b757a5e Mon Sep 17 00:00:00 2001 From: elmosa Date: Tue, 13 Oct 2026 12:24:08 +0200 Subject: [PATCH 24/25] edit changelog --- CHANGELOG.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d8fc460e8..04b17886e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,8 @@ All notable changes to this project will be documented in this file. ## Unreleased ### Added - +- #596 : Add `topetsc` method to `psydac.linalg.basic.ZeroOperator` +- #596 : Add `out` argument to `array_to_psydac` and `vec_topetsc` - #576 : Add module `psydac.feec.polar.conga_projections` with broken FEEC polar projections in 2D - #576 : Add 2D Poisson and TE Maxwell examples on polar mapped domains in `psydac/feec/polar/examples` - #576 : Add module `psydac.utilities.parallel_utils` for parallel execution and gathering variable-length arrays @@ -15,6 +16,8 @@ All notable changes to this project will be documented in this file. ### Fixed +- #596 : Allow export to VTK with a prescribed grid containing exact domain bounds +- #596 : Fix global geometric projectors with complex dtype - #576 : Fix bug in `TensorFemSpace.eval_field` caused by round-off at MPI subdomain boundaries - #576 : Require `sympde==0.19.3` which fixes a bug in the linearity checks - #579 : Require `h5py>=3.16` which installs correctly with `setuptools>=81.0` From 24bc46ec7f8ae0a38e5ae4ee880acaf0e1a5a3a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 13 Oct 2026 13:39:23 +0200 Subject: [PATCH 25/25] Apply small suggestion on CHANGELOG --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 04b17886e..48039f532 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -17,7 +17,7 @@ All notable changes to this project will be documented in this file. ### Fixed - #596 : Allow export to VTK with a prescribed grid containing exact domain bounds -- #596 : Fix global geometric projectors with complex dtype +- #596 : Fix global geometric projectors with complex `dtype` - #576 : Fix bug in `TensorFemSpace.eval_field` caused by round-off at MPI subdomain boundaries - #576 : Require `sympde==0.19.3` which fixes a bug in the linearity checks - #579 : Require `h5py>=3.16` which installs correctly with `setuptools>=81.0`