Merge branch 'devel' into test-pyccel-opt-flags

.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

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

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')

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

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

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