Small general fixes

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`

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:

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)
+    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 )
+    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)
+    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.5 )
+    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,12 +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.2 )
+    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)
@@ -181,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')
 
@@ -193,11 +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"])
+    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 * 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
@@ -253,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')
 
@@ -265,12 +291,20 @@ 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"])
+    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.3 )
-    f2  = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 )
+    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 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 )
 
     # Compute the projection
     u0 = P0(f1)
@@ -311,7 +345,7 @@ def test_derham_projector_2d_hcurl(ncells, degree, periodic, multiplicity, verbo
     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]])
@@ -324,6 +358,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<m doesn't work)
     multiplicity = [min(m, p) for p, m in zip(degree, multiplicity)]
 
@@ -333,7 +371,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 ) * 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,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")

psydac/fem/splines.py

@@ -185,7 +185,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 )
@@ -221,7 +221,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 )

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)
 

psydac/linalg/block.py

@@ -482,11 +482,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
 
 #===============================================================================
@@ -886,17 +886,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 ):
+    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):

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

psydac/linalg/tests/test_block.py

@@ -904,6 +904,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)
 
     # Initialize random number generator (set seed for reproducibility)
     rng = np.random.default_rng(seed=n1*n2*p1*p2)
@@ -917,10 +919,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] )
@@ -1343,6 +1351,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)
 
     # Initialize random number generator (set seed for reproducibility)
     rng = np.random.default_rng(seed=n1*n2*p1*p2)
@@ -1355,9 +1365,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.array_equal(x.toarray(), y.toarray())
+    assert np.array_equal(x.toarray(), z.toarray())
+    assert np.array_equal(v.array, u.array)
 
-    assert np.allclose( x.toarray() , v.toarray(), rtol=1e-12, atol=1e-12 )
 
 #=============================================================================== 
 @pytest.mark.parametrize( 'dtype', [float, complex] )