Set coordinates to real symbols (#170)

Domain coordinates in SymPDE are SymPy symbols. Here we set their
attribute `real=True` to avoid problems with code generation in Psydac
when working with complex expressions. Fixes #169.

sympde/expr/evaluation.py

@@ -567,7 +567,7 @@ def eval(cls, expr, domain):
                     if isinstance(domain, (NCube, NCubeInterior)):
                         bounds      = domain.min_coords if ext == -1 else domain.max_coords
                         J           = J.subs(coordinates[axis], bounds[axis])
-                    newcoords = [Symbol(c.name+"_plus") for c in coordinates]
+                    newcoords = [Symbol(c.name+"_plus", real=True) for c in coordinates]
                     subs      = list(zip(coordinates, newcoords))
                     J = J.subs(subs)
 

sympde/expr/tests/test_expr.py

@@ -35,6 +35,7 @@
 from sympde.expr.evaluation import TerminalExpr
 
 #==============================================================================
+
 def test_linear_expr_2d_1():
 
     domain = Domain('Omega', dim=2)
@@ -897,6 +898,30 @@ def test_terminal_expr_bilinear_3d_1():
     assert e2[0].expr          == dx(um)*dx(vm) + dy(um)*dy(vm) + dz(um)*dz(vm)
     assert e3[0].expr.factor() == (dx1(u)*dx1(v) + dx2(u)*dx2(v) + dx3(u)*dx3(v))*det
 
+#==============================================================================
+def test_terminal_expr_1d_1():
+    domain = Domain('Omega', dim=1)
+    x = domain.coordinates
+    assert TerminalExpr(grad(x), domain) == 1
+    assert TerminalExpr(grad(exp(1j*x) - x + 7), domain) == 1j*exp(1j*x) - 1
+
+#==============================================================================
+def test_terminal_expr_2d_1():
+    domain = Domain('Omega', dim=2)
+    x, y = domain.coordinates
+    assert TerminalExpr(grad(x), domain) == Matrix([[1], [0]])
+    assert TerminalExpr(grad(y), domain) == Matrix([[0], [1]])
+    assert TerminalExpr(grad(exp(1j*x)*sin(2*y)), domain) == Matrix([[1j*exp(1j*x)*sin(2*y)], [2*exp(1j*x)*cos(2*y)]])
+
+#==============================================================================
+def test_terminal_expr_3d_1():
+    domain = Domain('Omega', dim=3)
+    x, y, z = domain.coordinates
+    assert TerminalExpr(grad(x), domain) == Matrix([[1], [0], [0]])
+    assert TerminalExpr(grad(y), domain) == Matrix([[0], [1], [0]])
+    assert TerminalExpr(grad(z), domain) == Matrix([[0], [0], [1]])
+    assert TerminalExpr(grad(exp(1j*x)*sin(2*y)*z), domain) == Matrix([[1j*exp(1j*x)*sin(2*y)*z], [2*exp(1j*x)*cos(2*y)*z], [exp(1j*x)*sin(2*y)]])
+
 #==============================================================================
 @pytest.mark.skip(reason="New linearize() function does not accept 'LinearExpr' objects")
 def test_linearize_expr_2d_1():

sympde/exterior/tests/test_exterior.py

@@ -17,7 +17,7 @@
 #==============================================================================
 def test_exterior_1():
 
-    x, y, z = symbols('x y z')
+    x, y, z = symbols('x y z', real=True)
     a = Constant('a')
     n = Symbol('n')
 

sympde/topology/basic.py

@@ -32,7 +32,7 @@ def coordinates(self):
             else:
                 xyz = ['x', 'y', 'z'][:dim]
 
-            xyz = [Symbol(i) for i in xyz]
+            xyz = [Symbol(i, real=True) for i in xyz]
             self._coordinates = xyz
 
         if dim == 1:

sympde/topology/derivatives.py

@@ -104,7 +104,7 @@ def eval(cls, expr):
                 return S.Zero
 
             elif isinstance(expr, Expr):
-                x = Symbol(cls.coordinate)
+                x = Symbol(cls.coordinate, real=True)
                 if cls.logical:
                     M = expr.atoms(Mapping)
                     if len(M)>0:

sympde/topology/domain.py

@@ -839,7 +839,8 @@ def __new__(cls, name, dim, min_coords, max_coords):
             raise ValueError("Min coordinates must be smaller than max")
 
         coord_names = 'x1:{}'.format(dim + 1)
-        coordinates = symbols(coord_names)  # [MCP 24.07.2024] we should be able to pass real=True but requires additional testing as it breaks some calls to TerminalExpr
+
+        coordinates = symbols(coord_names, real=True)
 
         # Choose which type to use:
         #   a) if dim <= 3, use Line, Square or Cube;

sympde/topology/mapping.py

@@ -184,7 +184,7 @@ def __new__(cls, name, dim=None, **kwargs):
             return obj
 
         if coordinates is None:
-            _coordinates = [Symbol(name) for name in ['x', 'y', 'z'][:pdim]]
+            _coordinates = [Symbol(name, real=True) for name in ['x', 'y', 'z'][:pdim]]
         else:
             if not isinstance(coordinates, (list, tuple, Tuple)):
                 raise TypeError('> Expecting list, tuple, Tuple')
@@ -193,7 +193,7 @@ def __new__(cls, name, dim=None, **kwargs):
                 if not isinstance(a, (str, Symbol)):
                     raise TypeError('> Expecting str or Symbol')
 
-            _coordinates = [Symbol(u) for u in coordinates]
+            _coordinates = [Symbol(u, real=True) for u in coordinates]
 
         obj._name                = name
         obj._ldim                = ldim
@@ -203,16 +203,19 @@ def __new__(cls, name, dim=None, **kwargs):
         obj._is_minus            = None
         obj._is_plus             = None
 
-        lcoords = ['x1', 'x2', 'x3'][:ldim]
-        lcoords = [Symbol(i) for i in lcoords]
-        obj._logical_coordinates = Tuple(*lcoords)
+        lcoords_names = ['x1', 'x2', 'x3'][:ldim]
+        lcoords_symbols_real = [Symbol(i, real=True) for i in lcoords_names]
+        lcoords_symbols_real_dict = {key: symbol for key,symbol in zip(lcoords_names, lcoords_symbols_real)}
+        obj._logical_coordinates = Tuple(*lcoords_symbols_real) # coordinates must be symbols with real=True      
+        lcoords_general_symbols = [Symbol(i) for i in lcoords_names] #list of symbols without real=True       
+
         # ...
         if not( obj._expressions is None ):
-            coords = ['x', 'y', 'z'][:pdim]
+            coords_names = ['x', 'y', 'z'][:pdim]
 
             # ...
             args = []
-            for i in coords:
+            for i in coords_names:
                 x = obj._expressions[i]
                 x = sympify(x)
                 args.append(x)
@@ -225,13 +228,15 @@ def __new__(cls, name, dim=None, **kwargs):
                 x = sympify(i)
                 args = args.subs(x,0)
             # ...
-
-            constants        = list(set(args.free_symbols) - set(lcoords))
+            # get constants by subtracting coordinates from list of free symbols
+            constants        = list(set(args.free_symbols) - set(lcoords_general_symbols))
             constants_values = {a.name:Constant(a.name) for a in constants}
             # subs constants as Constant objects instead of Symbol
             constants_values.update( kwargs )
             d = {a:numpy_to_native_python(constants_values[a.name]) for a in constants}
             args = args.subs(d)
+            # subs coordinate symbols without real=True to symbols with real=True
+            args = args.subs(lcoords_symbols_real_dict)
 
             obj._expressions = args
             obj._constants   = tuple(a for a in constants if isinstance(constants_values[a.name], Symbol))
@@ -919,7 +924,7 @@ def eval(cls, expr, domain, **options):
             if has(expr, DiffOperator):
                 return cls( expr, domain, evaluate=False)
             else:
-                syms = symbols(ph_coords[:dim])
+                syms = symbols(ph_coords[:dim], real=True)
                 if isinstance(mapping, InterfaceMapping):
                     mapping = mapping.minus
                     # here we assume that the two mapped domains

sympde/topology/tests/test_gallery.py

@@ -4,7 +4,7 @@
 from sympy import Tuple
 from sympy import symbols
 
-x1, x2, x3 = symbols('x1, x2, x3')
+x1, x2, x3 = symbols('x1, x2, x3', real=True)
 
 from sympde.topology import Interval, ProductDomain, InteriorDomain, Domain
 from sympde.topology import Line, Square, Cube, NCubeInterior

sympde/topology/tests/test_logical_expr.py

@@ -585,7 +585,7 @@ def test_symbolic_expr_3d_1():
 def test_identity_mapping_2d_1():
     dim = 2
 
-    x1, x2 = symbols('x1, x2')
+    x1, x2 = symbols('x1, x2', real=True)
 
     M = IdentityMapping('M', dim=dim)
 
@@ -626,7 +626,7 @@ def test_identity_mapping_2d_2():
 def test_polar_mapping_2d_1():
     dim = 2
 
-    x1, x2 = symbols('x1, x2')
+    x1, x2 = symbols('x1, x2', real=True)
 
     constants = ['c1', 'c2', 'rmax', 'rmin']
     c1, c2, rmax, rmin = [Constant(i) for i in constants]
@@ -656,7 +656,7 @@ def test_polar_mapping_2d_1():
 def test_target_mapping_2d_1():
     dim = 2
 
-    x1, x2 = symbols('x1, x2')
+    x1, x2 = symbols('x1, x2', real=True)
 
     constants = ['c1', 'c2', 'D', 'k']
     c1, c2, D, k = [Constant(i) for i in constants]
@@ -688,7 +688,7 @@ def test_target_mapping_2d_1():
 def test_czarny_mapping_2d_1():
     dim = 2
 
-    x1, x2 = symbols('x1, x2')
+    x1, x2 = symbols('x1, x2', real=True)
 
     constants = ['c2', 'eps', 'b']
     c2, eps, b = [Constant(i) for i in constants]
@@ -725,7 +725,7 @@ def test_czarny_mapping_2d_1():
 def test_collela_mapping_2d_1():
     dim = 2
 
-    x1, x2 = symbols('x1, x2')
+    x1, x2 = symbols('x1, x2', real=True)
 
     constants = ['eps', 'k1', 'k2']
     eps, k1, k2 = [Constant(i) for i in constants]
@@ -763,7 +763,7 @@ def test_collela_mapping_2d_1():
 def test_torus_mapping_3d_1():
     dim = 3
 
-    x1, x2, x3 = symbols('x1, x2, x3')
+    x1, x2, x3 = symbols('x1, x2, x3', real=True)
     R0 = Constant('R0')
 
     M = TorusMapping('M', dim=dim)
@@ -812,7 +812,7 @@ def test_torus_mapping_3d_1():
 def test_twisted_target_mapping_3d_1():
     dim = 3
 
-    x1, x2, x3 = symbols('x1, x2, x3')
+    x1, x2, x3 = symbols('x1, x2, x3', real=True)
 
     constants = ['c1', 'c2', 'c3', 'D', 'k']
     c1, c2, c3, D, k = [Constant(i) for i in constants]