Code generation fails when taking real or imaginary part of symbolic expression

[e-moral-sanchez]

The generated code does not work well when using Sympy functions re and im . There are two issues:

1. In the generated code Pyccel does not understand the keywords re and im.
2. The printer prints the real and imaginary of coordinates, which should be always real.

Minimal example that crashes:

from sympde.expr                    import LinearForm, TerminalExpr, integral
from sympde.calculus                import dot
from sympde.topology                import Cube, Derham, element_of
from psydac.api.settings            import PSYDAC_BACKENDS
from psydac.api.discretization      import discretize
import sympy as sy

domain = Cube('domain')
derham = Derham(domain)

x,y,z = domain.coordinates

domain_h = discretize(domain, ncells=[5,5,5])
derham_h = discretize(derham, domain_h, degree=[1,1,1])
derham.V0.codomain_type='complex'
derham.V1.codomain_type='complex'
derham.V2.codomain_type='complex'
derham.V3.codomain_type='complex'
print('domain and space discretized')

f = TerminalExpr(sy.Matrix([0,0, y]), domain)
v = element_of(derham.V1, name='v')
b = LinearForm(v, integral(domain, dot(sy.re(f),v)))

backend = PSYDAC_BACKENDS['pyccel-gcc']
b_h = discretize(b, domain_h, derham_h.V1, backend=backend).assemble()

Error:

ERROR at annotation (semantic) stage
pyccel:
 |fatal [semantic]: dependencies_hcvpx9db_x61q777yyolm.py [46,73]| Undefined function (re)

---------------------------------------------------------------------------
PyccelSemanticError                       Traceback (most recent call last)
<ipython-input-24-844cd901adb8> in <module>
     25 
     26 backend = PSYDAC_BACKENDS['pyccel-gcc']
---> 27 b_h = discretize(b, domain_h, derham_h.V1, backend=backend).assemble()
     28 
     29 print('done!')

~/psydac/psydac/api/discretization.py in discretize(a, *args, **kwargs)
    585 
    586     elif isinstance(a, sym_LinearForm):
--> 587         return DiscreteLinearForm(a, kernel_expr, *args, **kwargs)
    588 
    589     elif isinstance(a, sym_Functional):

~/psydac/psydac/api/fem.py in __init__(self, expr, kernel_expr, domain_h, space, nquads, vector, update_ghost_regions, backend, symbolic_mapping)
   1084         # BasicDiscrete generates the assembly code and sets the following attributes that are used afterwards:
   1085         # self._func, self._free_args, self._max_nderiv and self._backend
-> 1086         BasicDiscrete.__init__(self, expr, kernel_expr, comm=comm, root=0, discrete_space=discrete_space,
   1087                               nquads=nquads, is_rational_mapping=is_rational_mapping, mapping=symbolic_mapping,
   1088                               mapping_space=mapping_space, num_threads=self._num_threads, backend=backend)

~/psydac/psydac/api/basic.py in __init__(self, expr, kernel_expr, folder, comm, root, discrete_space, nquads, is_rational_mapping, mapping, mapping_space, num_threads, backend)
    296                        mapping_space=None, num_threads=None, backend=None):
    297 
--> 298         BasicCodeGen.__init__(self, expr, folder=folder, comm=comm, root=root, discrete_space=discrete_space,
    299                        kernel_expr=kernel_expr, nquads=nquads, is_rational_mapping=is_rational_mapping,
    300                        mapping=mapping, mapping_space=mapping_space, num_threads=num_threads, backend=backend)

~/psydac/psydac/api/basic.py in __init__(self, expr, folder, comm, root, discrete_space, kernel_expr, nquads, is_rational_mapping, mapping, mapping_space, num_threads, backend)
    148         if comm is not None and comm.size>1: comm.Barrier()
    149         # compile code
--> 150         self._compile()
    151 
    152     @property

~/psydac/psydac/api/basic.py in _compile(self)
    280 
    281         if self.backend['name'] == 'pyccel':
--> 282             package = self._compile_pyccel(package)
    283         elif self.backend['name'] == 'pythran':
    284             package = self._compile_pythran(package)

~/psydac/psydac/api/basic.py in _compile_pyccel(self, mod, verbose)
    261 
    262         from pyccel.epyccel import epyccel
--> 263         fmod = epyccel(mod,
    264                        accelerators = accelerators,
    265                        compiler    = compiler,

~/vpsydac4/lib/python3.10/site-packages/pyccel/epyccel.py in epyccel(python_function_or_module, **kwargs)
    387             mod, fun = epyccel_seq( python_function_or_module, **kwargs )
    388         except PyccelError as e:
--> 389             raise type(e)(str(e)) from None
    390 
    391     # Return Fortran function (if any), otherwise module

PyccelSemanticError: Semantic step failed

[yguclu]

Can you please share the generated Python code in the __psydac__ directory?

[e-moral-sanchez]

The file dependencies_hcvpx9db.py in __psydac__ directory is:

  # Not supported in Python:
  # re
def assemble_vector_hcvpx9db(global_test_basis_v_0_1 : "float64[:,:,:,:]", global_test_basis_v_0_2 : "float64[:,:,:,:]", global_test_basis_v_0_3 : "float64[:,:,:,:]", global_test_basis_v_1_1 : "float64[:,:,:,:]", global_test_basis_v_1_2 : "float64[:,:,:,:]", global_test_basis_v_1_3 : "float64[:,:,:,:]", global_test_basis_v_2_1 : "float64[:,:,:,:]", global_test_basis_v_2_2 : "float64[:,:,:,:]", global_test_basis_v_2_3 : "float64[:,:,:,:]", global_span_v_0_1 : "int64[:]", global_span_v_0_2 : "int64[:]", global_span_v_0_3 : "int64[:]", global_span_v_1_1 : "int64[:]", global_span_v_1_2 : "int64[:]", global_span_v_1_3 : "int64[:]", global_span_v_2_1 : "int64[:]", global_span_v_2_2 : "int64[:]", global_span_v_2_3 : "int64[:]", global_x1 : "float64[:,:]", global_x2 : "float64[:,:]", global_x3 : "float64[:,:]", test_v_0_p1 : "int64", test_v_0_p2 : "int64", test_v_0_p3 : "int64", test_v_1_p1 : "int64", test_v_1_p2 : "int64", test_v_1_p3 : "int64", test_v_2_p1 : "int64", test_v_2_p2 : "int64", test_v_2_p3 : "int64", n_element_1 : "int64", n_element_2 : "int64", n_element_3 : "int64", k1 : "int64", k2 : "int64", k3 : "int64", pad1 : "int64", pad2 : "int64", pad3 : "int64", g_vec_v_2_hcvpx9db : "complex[:,:,:]"):

    from numpy import array, zeros, zeros_like, floor
    local_x1 = zeros_like(global_x1[0,:])
    local_x2 = zeros_like(global_x2[0,:])
    local_x3 = zeros_like(global_x3[0,:])

    l_vec_v_2_hcvpx9db = zeros((2, 2, 1), dtype='complex')
    for i_element_1 in range(0, n_element_1, 1):
        local_x1[:] = global_x1[i_element_1,:]
        span_v_0_1 = global_span_v_0_1[i_element_1]
        span_v_1_1 = global_span_v_1_1[i_element_1]
        span_v_2_1 = global_span_v_2_1[i_element_1]
        for i_element_2 in range(0, n_element_2, 1):
            local_x2[:] = global_x2[i_element_2,:]
            span_v_0_2 = global_span_v_0_2[i_element_2]
            span_v_1_2 = global_span_v_1_2[i_element_2]
            span_v_2_2 = global_span_v_2_2[i_element_2]
            for i_element_3 in range(0, n_element_3, 1):
                local_x3[:] = global_x3[i_element_3,:]
                span_v_0_3 = global_span_v_0_3[i_element_3]
                span_v_1_3 = global_span_v_1_3[i_element_3]
                span_v_2_3 = global_span_v_2_3[i_element_3]
                for i_basis_1 in range(0, 1 + test_v_2_p1, 1):
                    for i_basis_2 in range(0, 1 + test_v_2_p2, 1):
                        for i_basis_3 in range(0, 1 + test_v_2_p3, 1):
                            contribution_v_2_hcvpx9db = 0j
                            for i_quad_1 in range(0, 2, 1):
                                x1 = local_x1[i_quad_1]
                                v_2_1 = global_test_basis_v_2_1[i_element_1,i_basis_1,0,i_quad_1]
                                v_2_1_x1 = global_test_basis_v_2_1[i_element_1,i_basis_1,1,i_quad_1]
                                for i_quad_2 in range(0, 2, 1):
                                    x2 = local_x2[i_quad_2]
                                    v_2_2 = global_test_basis_v_2_2[i_element_2,i_basis_2,0,i_quad_2]
                                    v_2_2_x2 = global_test_basis_v_2_2[i_element_2,i_basis_2,1,i_quad_2]
                                    for i_quad_3 in range(0, 2, 1):
                                        x3 = local_x3[i_quad_3]
                                        v_2_3 = global_test_basis_v_2_3[i_element_3,i_basis_3,0,i_quad_3]
                                        v_2_3_x3 = global_test_basis_v_2_3[i_element_3,i_basis_3,1,i_quad_3]
                                        v_2 = v_2_1*v_2_2*v_2_3
                                        v_2_x3 = v_2_1*v_2_2*v_2_3_x3
                                        v_2_x2 = v_2_1*v_2_2_x2*v_2_3
                                        v_2_x1 = v_2_1_x1*v_2_2*v_2_3
                                        contribution_v_2_hcvpx9db += v_2*re(x2)

                            l_vec_v_2_hcvpx9db[i_basis_1,i_basis_2,i_basis_3] = contribution_v_2_hcvpx9db

                g_vec_v_2_hcvpx9db[pad1 + span_v_2_1 - test_v_2_p1:1 + pad1 + span_v_2_1,pad2 + span_v_2_2 - test_v_2_p2:1 + pad2 + span_v_2_2,pad3 + span_v_2_3 - test_v_2_p3:1 + pad3 + span_v_2_3] += l_vec_v_2_hcvpx9db[:,:,:]

    return

[yguclu]

Indeed, it appears that Psydac does not know how to map SymPy's re function to the method .real which is native for Python scalar numeric types, and also available in NumPy arrays.

[yguclu]

In other words, re(x2) should be replaced with x2.real in the code snippet above.