diff --git a/.github/workflows/continuous-integration.yml b/.github/workflows/continuous-integration.yml
index ee8ebcfab..7eedd019e 100644
--- a/.github/workflows/continuous-integration.yml
+++ b/.github/workflows/continuous-integration.yml
@@ -16,14 +16,17 @@ jobs:
       fail-fast: false
       matrix:
         os: [ ubuntu-latest ]
-        python-version: [ 3.9, '3.10', '3.11' ]
+        python-version: [ 3.9, '3.10', '3.11', '3.12' ]
         isMerge:
           - ${{ github.event_name == 'push' && github.ref == 'refs/heads/devel' }}
         exclude:
           - { isMerge: false, python-version: '3.10' }
+          - { isMerge: false, python-version: '3.11' }
         include:
           - os: macos-latest
             python-version: '3.10'
+          - os: macos-latest
+            python-version: '3.11'
 
     name: ${{ matrix.os }} / Python ${{ matrix.python-version }}
 
@@ -95,7 +98,7 @@ jobs:
 
       - name: Download a specific release of PETSc
         run: |
-          git clone --depth 1 --branch v3.20.5 https://gitlab.com/petsc/petsc.git
+          git clone --depth 1 --branch v3.21.4 https://gitlab.com/petsc/petsc.git
 
       - name: Install PETSc with complex support, and test it
         working-directory: ./petsc
@@ -115,7 +118,8 @@ jobs:
 
       - name: Install Python dependencies
         run: |
-          export CC="mpicc" HDF5_MPI="ON"
+          export CC="mpicc"
+          export HDF5_MPI="ON"
           python -m pip install -r requirements.txt
           python -m pip install -r requirements_extra.txt --no-build-isolation
           python -m pip list
diff --git a/README.md b/README.md
index f674cd290..c35348982 100644
--- a/README.md
+++ b/README.md
@@ -131,7 +131,7 @@ Although Psydac provides several iterative linear solvers which work with our na
 In order to use these additional feature, PETSc and petsc4py must be installed as follows.
 First, we download the latest release of PETSc from its [official Git repository](https://gitlab.com/petsc/petsc):
 ```sh
-git clone --depth 1 --branch v3.20.5 https://gitlab.com/petsc/petsc.git
+git clone --depth 1 --branch v3.21.4 https://gitlab.com/petsc/petsc.git
 ```
 Next, we specify a configuration for complex numbers, and install PETSc in a local directory:
 ```sh
diff --git a/examples/maxwell_2d_multi_patch.py b/examples/maxwell_2d_multi_patch.py
index 915bd214a..b731d18be 100644
--- a/examples/maxwell_2d_multi_patch.py
+++ b/examples/maxwell_2d_multi_patch.py
@@ -129,11 +129,12 @@ def run_maxwell_2d(uex, f, alpha, domain, ncells, degree, k=None, kappa=None, co
     N = 20
 
     mappings = OrderedDict([(P.logical_domain, P.mapping) for P in domain.interior])
-    mappings_list = [m.get_callable_mapping() for m in mappings.values()]
+    mappings_list = [m for m in mappings.values()]
+    call_mappings_list = [m.get_callable_mapping() for m in mappings_list]
 
     Eex_x   = lambdify(domain.coordinates, Eex[0])
     Eex_y   = lambdify(domain.coordinates, Eex[1])
-    Eex_log = [pull_2d_hcurl([Eex_x, Eex_y], f) for f in mappings_list]
+    Eex_log = [pull_2d_hcurl([Eex_x, Eex_y], f) for f in call_mappings_list]
 
     etas, xx, yy    = get_plotting_grid(mappings, N=20)
     grid_vals_hcurl = lambda v: get_grid_vals(v, etas, mappings_list, space_kind='hcurl')
diff --git a/psydac/api/tests/test_2d_complex.py b/psydac/api/tests/test_2d_complex.py
index 57e6a37db..9caf559cc 100644
--- a/psydac/api/tests/test_2d_complex.py
+++ b/psydac/api/tests/test_2d_complex.py
@@ -1,11 +1,14 @@
 # -*- coding: UTF-8 -*-
 import os
 
+from collections import OrderedDict
+
 from mpi4py import MPI
 import pytest
 import numpy as np
 from sympy import pi, cos, sin, symbols, conjugate, exp
 from sympy import Tuple, Matrix
+from sympy import lambdify
 
 from sympde.calculus import grad, dot, cross, curl
 from sympde.calculus import minus, plus
@@ -244,7 +247,7 @@ def run_helmholtz_2d(solution, kappa, e_w_0, dx_e_w_0, domain, ncells=None, degr
     l2_error = l2norm_h.assemble(u=uh)
     h1_error = h1norm_h.assemble(u=uh)
 
-    return l2_error, h1_error
+    return l2_error, h1_error, uh
 
 #==============================================================================
 def run_maxwell_2d(uex, f, alpha, domain, *, ncells=None, degree=None, filename=None, comm=None):
@@ -400,7 +403,7 @@ def test_complex_poisson_2d_multipatch_mapping():
     assert abs(h1_error - expected_h1_error) < 1e-7
 
 
-def test_complex_helmholtz_2d():
+def test_complex_helmholtz_2d(plot_sol=False):
     # This test solves the homogeneous Helmhotz equation with impedance BC. 
     # In particular, we impose an incoming wave from the left and absorbing boundary conditions at the right.
     # Along y, periodic boundary conditions are considered.
@@ -412,11 +415,50 @@ def test_complex_helmholtz_2d():
     e_w_0 = sin(kappa * y) # value of incoming wave at x=0, forall y
     dx_e_w_0 = 1j*kappa*sin(kappa * y) # derivative wrt. x of incoming wave at x=0, forall y
 
-    l2_error, h1_error = run_helmholtz_2d(solution, kappa, e_w_0, dx_e_w_0, domain, ncells=[2**2,2**2], degree=[2,2])
+    l2_error, h1_error, uh = run_helmholtz_2d(solution, kappa, e_w_0, dx_e_w_0, domain, ncells=[2**2,2**2], degree=[2,2])
 
     expected_l2_error = 0.01540947560953227
     expected_h1_error = 0.19040207344639598
 
+    print(f'errors: l2 = {l2_error}, h1 = {h1_error}')
+    print('expected errors: l2 = {}, h1 = {}'.format(expected_l2_error, expected_h1_error))
+    
+    if plot_sol:
+        from psydac.feec.multipatch.plotting_utilities import get_plotting_grid, get_grid_vals
+        from psydac.feec.multipatch.plotting_utilities import get_patch_knots_gridlines, my_small_plot
+        from psydac.feec.pull_push                     import pull_2d_h1
+        
+        Id_mapping = IdentityMapping('M', 2)
+        # print(f'domain.interior = {domain.interior}')
+        # domain_interior = [domain]
+        # print(f'domain.logical_domain = {domain.logical_domain}')
+        mappings = OrderedDict([(domain, Id_mapping)])
+        mappings_list = [m for m in mappings.values()]
+        call_mappings_list = [m.get_callable_mapping() for m in mappings_list]
+
+        uh = [uh]  # single-patch cast as multi-patch solution 
+
+        u   = lambdify(domain.coordinates, solution)
+        u_log = [pull_2d_h1(u, f) for f in call_mappings_list]
+
+        etas, xx, yy         = get_plotting_grid(mappings, N=20)
+        grid_vals_h1         = lambda v: get_grid_vals(v, etas, mappings_list, space_kind='h1')
+
+        uh_vals = grid_vals_h1(uh)
+        u_vals  = grid_vals_h1(u_log)
+
+        u_err   = [(u1 - u2) for u1, u2 in zip(u_vals, uh_vals)]
+    
+        my_small_plot(
+            title=r'approximation of solution $u$',
+            vals=[u_vals, uh_vals, u_err],
+            titles=[r'$u^{ex}(x,y)$', r'$u^h(x,y)$', r'$|(u^{ex}-u^h)(x,y)|$'],
+            xx=xx,
+            yy=yy,
+            gridlines_x1=None,
+            gridlines_x2=None,
+        )
+
     assert( abs(l2_error - expected_l2_error) < 1.e-7)
     assert( abs(h1_error - expected_h1_error) < 1.e-7)
 
@@ -487,59 +529,62 @@ def teardown_function():
 
 if __name__ == '__main__':
 
-    from collections import OrderedDict
-
-    from sympy       import lambdify
-
-    from psydac.feec.multipatch.plotting_utilities import get_plotting_grid, get_grid_vals
-    from psydac.feec.multipatch.plotting_utilities import get_patch_knots_gridlines, my_small_plot
-    from psydac.api.tests.build_domain             import build_pretzel
-    from psydac.feec.pull_push                     import pull_2d_hcurl
+    check = 'complex_helmholtz_2d'
 
-    domain = build_pretzel()
-    x,y    = domain.coordinates
+    if check == 'complex_helmholtz_2d':
+        test_complex_helmholtz_2d(plot_sol=True)
 
-    omega = 1.5
-    alpha = -omega**2
-    Eex   = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))
-    f     = Tuple(alpha*sin(pi*y) - pi**2*sin(pi*y)*cos(pi*x) + pi**2*sin(pi*y),
-                  alpha*sin(pi*x)*cos(pi*y) + pi**2*sin(pi*x)*cos(pi*y))
+    else:
 
-    l2_error, Eh = run_maxwell_2d(Eex, f, alpha, domain, ncells=[2**2, 2**2], degree=[2,2])
-
-    mappings = OrderedDict([(P.logical_domain, P.mapping) for P in domain.interior])
-    mappings_list = list(mappings.values())
-    mappings_list = [mapping.get_callable_mapping() for mapping in mappings_list]
-
-    Eex_x   = lambdify(domain.coordinates, Eex[0])
-    Eex_y   = lambdify(domain.coordinates, Eex[1])
-    Eex_log = [pull_2d_hcurl([Eex_x,Eex_y], f) for f in mappings_list]
-
-    etas, xx, yy         = get_plotting_grid(mappings, N=20)
-    grid_vals_hcurl      = lambda v: get_grid_vals(v, etas, mappings_list, space_kind='hcurl')
-
-    Eh_x_vals, Eh_y_vals = grid_vals_hcurl(Eh)
-    E_x_vals, E_y_vals   = grid_vals_hcurl(Eex_log)
-
-    E_x_err              = [(u1 - u2) for u1, u2 in zip(E_x_vals, Eh_x_vals)]
-    E_y_err              = [(u1 - u2) for u1, u2 in zip(E_y_vals, Eh_y_vals)]
-
-    my_small_plot(
-        title=r'approximation of solution $u$, $x$ component',
-        vals=[E_x_vals, Eh_x_vals, E_x_err],
-        titles=[r'$u^{ex}_x(x,y)$', r'$u^h_x(x,y)$', r'$|(u^{ex}-u^h)_x(x,y)|$'],
-        xx=xx,
-        yy=yy,
-        gridlines_x1=None,
-        gridlines_x2=None,
-    )
-
-    my_small_plot(
-        title=r'approximation of solution $u$, $y$ component',
-        vals=[E_y_vals, Eh_y_vals, E_y_err],
-        titles=[r'$u^{ex}_y(x,y)$', r'$u^h_y(x,y)$', r'$|(u^{ex}-u^h)_y(x,y)|$'],
-        xx=xx,
-        yy=yy,
-        gridlines_x1=None,
-        gridlines_x2=None,
-    )
+        from psydac.feec.multipatch.plotting_utilities import get_plotting_grid, get_grid_vals
+        from psydac.feec.multipatch.plotting_utilities import get_patch_knots_gridlines, my_small_plot
+        from psydac.api.tests.build_domain             import build_pretzel
+        from psydac.feec.pull_push                     import pull_2d_hcurl
+        
+        domain = build_pretzel()
+        x,y    = domain.coordinates
+
+        omega = 1.5
+        alpha = -omega**2
+        Eex   = Tuple(sin(pi*y), sin(pi*x)*cos(pi*y))
+        f     = Tuple(alpha*sin(pi*y) - pi**2*sin(pi*y)*cos(pi*x) + pi**2*sin(pi*y),
+                    alpha*sin(pi*x)*cos(pi*y) + pi**2*sin(pi*x)*cos(pi*y))
+
+        l2_error, Eh = run_maxwell_2d(Eex, f, alpha, domain, ncells=[2**2, 2**2], degree=[2,2])
+
+        mappings = OrderedDict([(P.logical_domain, P.mapping) for P in domain.interior])
+        mappings_list = list(mappings.values())
+        call_mappings_list = [m.get_callable_mapping() for m in mappings_list]
+
+        Eex_x   = lambdify(domain.coordinates, Eex[0])
+        Eex_y   = lambdify(domain.coordinates, Eex[1])
+        Eex_log = [pull_2d_hcurl([Eex_x,Eex_y], f) for f in call_mappings_list]
+
+        etas, xx, yy         = get_plotting_grid(mappings, N=20)
+        grid_vals_hcurl      = lambda v: get_grid_vals(v, etas, mappings_list, space_kind='hcurl')
+
+        Eh_x_vals, Eh_y_vals = grid_vals_hcurl(Eh)
+        E_x_vals, E_y_vals   = grid_vals_hcurl(Eex_log)
+
+        E_x_err              = [(u1 - u2) for u1, u2 in zip(E_x_vals, Eh_x_vals)]
+        E_y_err              = [(u1 - u2) for u1, u2 in zip(E_y_vals, Eh_y_vals)]
+
+        my_small_plot(
+            title=r'approximation of solution $u$, $x$ component',
+            vals=[E_x_vals, Eh_x_vals, E_x_err],
+            titles=[r'$u^{ex}_x(x,y)$', r'$u^h_x(x,y)$', r'$|(u^{ex}-u^h)_x(x,y)|$'],
+            xx=xx,
+            yy=yy,
+            gridlines_x1=None,
+            gridlines_x2=None,
+        )
+
+        my_small_plot(
+            title=r'approximation of solution $u$, $y$ component',
+            vals=[E_y_vals, Eh_y_vals, E_y_err],
+            titles=[r'$u^{ex}_y(x,y)$', r'$u^h_y(x,y)$', r'$|(u^{ex}-u^h)_y(x,y)|$'],
+            xx=xx,
+            yy=yy,
+            gridlines_x1=None,
+            gridlines_x2=None,
+        )
diff --git a/psydac/api/tests/test_2d_navier_stokes.py b/psydac/api/tests/test_2d_navier_stokes.py
index 867e6cd8d..a0bea91d4 100644
--- a/psydac/api/tests/test_2d_navier_stokes.py
+++ b/psydac/api/tests/test_2d_navier_stokes.py
@@ -362,15 +362,15 @@ def test_st_navier_stokes_2d():
 
     a = TerminalExpr(-mu*laplace(ue), domain)
     b = TerminalExpr(    grad(ue), domain)
-    c = TerminalExpr(    grad(pe), domain)
+    c = TerminalExpr(    grad(pe), domain)    
     f = (a + b.T*ue + c).simplify()
-
+    
     fx = -mu*(ux.diff(x, 2) + ux.diff(y, 2)) + ux*ux.diff(x) + uy*ux.diff(y) + pe.diff(x)
-    fy = -mu*(uy.diff(x, 2) - uy.diff(y, 2)) + ux*uy.diff(x) + uy*uy.diff(y) + pe.diff(y)
+    fy = -mu*(uy.diff(x, 2) + uy.diff(y, 2)) + ux*uy.diff(x) + uy*uy.diff(y) + pe.diff(y)
 
     assert (f[0]-fx).simplify() == 0
     assert (f[1]-fy).simplify() == 0
-
+    
     f  = Tuple(fx, fy)
     # ...
 
@@ -420,7 +420,7 @@ def test_st_navier_stokes_2d_parallel():
     f = (a + b.T*ue + c).simplify()
 
     fx = -mu*(ux.diff(x, 2) + ux.diff(y, 2)) + ux*ux.diff(x) + uy*ux.diff(y) + pe.diff(x)
-    fy = -mu*(uy.diff(x, 2) - uy.diff(y, 2)) + ux*uy.diff(x) + uy*uy.diff(y) + pe.diff(y)
+    fy = -mu*(uy.diff(x, 2) + uy.diff(y, 2)) + ux*uy.diff(x) + uy*uy.diff(y) + pe.diff(y)
 
     assert (f[0]-fx).simplify() == 0
     assert (f[1]-fy).simplify() == 0
@@ -450,20 +450,30 @@ def teardown_function():
 
 #------------------------------------------------------------------------------
 if __name__ == '__main__':
-    import matplotlib.pyplot as plt
-    from matplotlib import animation
-    from time       import time
 
-    Tf       = 3.
-    dt_h     = 0.05
-    nt       = Tf//dt_h
-    filename = os.path.join(mesh_dir, 'bent_pipe.h5')
+    verify = 'st_navier_stokes_2d'
+
+    if verify == 'st_navier_stokes_2d':
+        print('Running test_st_navier_stokes_2d_parallel()')
+        test_st_navier_stokes_2d_parallel()
+        # test_st_navier_stokes_2d()
+    
+    else:
+    
+        import matplotlib.pyplot as plt
+        from matplotlib import animation
+        from time       import time
+
+        Tf       = 3.
+        dt_h     = 0.05
+        nt       = Tf//dt_h
+        filename = os.path.join(mesh_dir, 'bent_pipe.h5')
 
-    solutions, p_h, domain, domain_h = run_time_dependent_navier_stokes_2d(filename, dt_h=dt_h, nt=nt, scipy=False)
+        solutions, p_h, domain, domain_h = run_time_dependent_navier_stokes_2d(filename, dt_h=dt_h, nt=nt, scipy=False)
 
-    domain = domain.logical_domain
-    mapping = domain_h.mappings['patch_0']
+        domain = domain.logical_domain
+        mapping = domain_h.mappings['patch_0']
 
-    anim = animate_field(solutions, domain, mapping, res=(150,150), progress=True)
-    anim.save('animated_fields_{}_{}.mp4'.format(str(Tf).replace('.','_'), str(dt_h).replace('.','_')), writer=animation.FFMpegWriter(fps=60))
+        anim = animate_field(solutions, domain, mapping, res=(150,150), progress=True)
+        anim.save('animated_fields_{}_{}.mp4'.format(str(Tf).replace('.','_'), str(dt_h).replace('.','_')), writer=animation.FFMpegWriter(fps=60))
 
diff --git a/psydac/api/tests/test_api_2d_fields.py b/psydac/api/tests/test_api_2d_fields.py
index 9946af354..89f822ef6 100644
--- a/psydac/api/tests/test_api_2d_fields.py
+++ b/psydac/api/tests/test_api_2d_fields.py
@@ -153,7 +153,7 @@ def run_boundary_field_test(domain, boundary, f, ncells):
     a_grad_F_h = discretize(a_grad_F, domain_h, [Vh, Vh])
     a_grad_f_h = discretize(a_grad_f, domain_h, [Vh, Vh])
     a_grad_F_x = a_grad_F_h.assemble(F=fh)
-    a_grad_f_x = a_grad_f_h.assemble()    
+    a_grad_f_x = a_grad_f_h.assemble()
 
     lF_h = discretize(lF, domain_h,  Vh)
     lf_h = discretize(lf, domain_h,  Vh)
@@ -288,8 +288,8 @@ def test_field_identity_1():
     expected_error_1 =  4.77987181085604e-12
     expected_error_2 =  1.196388887893425e-07
 
-    assert( abs(error_1 - expected_error_1) < 1.e-7)
-    assert( abs(error_2 - expected_error_2) < 1.e-7)
+    assert abs(error_1 - expected_error_1) < 1.e-7
+    assert abs(error_2 - expected_error_2) < 1.e-7
 
 #------------------------------------------------------------------------------
 def test_field_identity_2():
@@ -302,8 +302,8 @@ def test_field_identity_2():
     expected_error_1 =  5.428295909559039e-11
     expected_error_2 =  2.9890068935570224e-11
 
-    assert( abs(error_1 - expected_error_1) < 1.e-10)
-    assert( abs(error_2 - expected_error_2) < 1.e-10)
+    assert abs(error_1 - expected_error_1) < 1.e-10
+    assert abs(error_2 - expected_error_2) < 1.e-10
 
 #------------------------------------------------------------------------------
 def test_field_collela():
@@ -316,8 +316,8 @@ def test_field_collela():
     expected_error_1 =  1.9180860719170134e-10
     expected_error_2 =  0.00010748308338081464
 
-    assert( abs(error_1 - expected_error_1) < 1.e-7)
-    assert( abs(error_2 - expected_error_2) < 1.e-7)
+    assert abs(error_1 - expected_error_1) < 1.e-7
+    assert abs(error_2 - expected_error_2) < 1.e-7
 
 #------------------------------------------------------------------------------
 def test_field_quarter_annulus():
@@ -332,8 +332,8 @@ def test_field_quarter_annulus():
     expected_error_1 =  1.1146377538410329e-10
     expected_error_2 =  9.18920469410037e-08
 
-    assert( abs(error_1 - expected_error_1) < 1.e-7)
-    assert( abs(error_2 - expected_error_2) < 1.e-7)
+    assert abs(error_1 - expected_error_1) < 1.e-7
+    assert abs(error_2 - expected_error_2) < 1.e-7
 
 #==============================================================================
 def test_nonlinear_poisson_circle():
@@ -343,7 +343,7 @@ def test_nonlinear_poisson_circle():
 
     expected_l2_error =  0.004026218710733066
 
-    assert( abs(l2_error - expected_l2_error) < 1.e-7)
+    assert abs(l2_error - expected_l2_error) < 1.e-7
 
 #==============================================================================
 # CLEAN UP SYMPY NAMESPACE
diff --git a/psydac/api/tests/test_assembly.py b/psydac/api/tests/test_assembly.py
index 853935199..44979df1b 100644
--- a/psydac/api/tests/test_assembly.py
+++ b/psydac/api/tests/test_assembly.py
@@ -1,8 +1,7 @@
 import pytest
 import numpy as np
 from mpi4py import MPI
-from sympy import pi, sin, cos, tan, atan, atan2, exp, sinh, cosh, tanh, atanh, Tuple, I
-
+from sympy import pi, sin, cos, tan, atan, atan2, exp, sinh, cosh, tanh, atanh, Tuple, I, sqrt
 
 from sympde.topology import Line, Square
 from sympde.topology import ScalarFunctionSpace, VectorFunctionSpace
@@ -219,10 +218,10 @@ def test_Norm_complex(backend):
     v = element_of(V, name='v')
     c = Constant(name='c', complex=True)
 
-    res = (1.+1.j)/np.sqrt(2)
+    res = (1.+1j)/np.sqrt(2)
 
     # We try to put complex as a sympy object in the expression
-    g1  = (1.+I)/np.sqrt(2)
+    g1  = (1.+I)/sqrt(2)
 
     # We try to put complex as a python scalar in the expression
     g2  = res
@@ -532,6 +531,8 @@ def test_assembly_no_synchr_args(backend):
 
 #==============================================================================
 if __name__ == '__main__':
+    test_Norm_complex(None)
+    exit()
     test_field_and_constant(None)
     test_multiple_fields(None)
     test_math_imports(None)
diff --git a/psydac/fem/tensor.py b/psydac/fem/tensor.py
index bc72578de..8832ad1a8 100644
--- a/psydac/fem/tensor.py
+++ b/psydac/fem/tensor.py
@@ -693,6 +693,10 @@ def integral(self, f, *, nquads=None):
             mpi_comm = self.vector_space.cart.comm
             c = mpi_comm.allreduce(c)
 
+        # convert to native python type if numpy to avoid errors with sympify
+        if isinstance(c, np.generic):
+            c = c.item()
+        
         return c
 
     #--------------------------------------------------------------------------
diff --git a/pyproject.toml b/pyproject.toml
index a797463b5..032ea1dda 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -7,7 +7,7 @@ name            = "psydac"
 version         = "0.1"
 description     = "Python package for isogeometric analysis (IGA)"
 readme          = "README.md"
-requires-python = ">= 3.9, < 3.12"
+requires-python = ">= 3.9, < 3.13"
 license         = {file = "LICENSE"}
 authors         = [
     {name = "Psydac development team", email = "psydac@googlegroups.com"}
@@ -33,7 +33,7 @@ dependencies = [
     'pyevtk',
 
     # Our packages from PyPi
-    'sympde == 0.18.3',
+    'sympde == 0.19.0',
     'pyccel >= 1.11.2',
     'gelato == 0.12',