From 9f115bf9e4a65f4e82dd870114ac8a338ffbca49 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 26 Aug 2019 18:00:58 +0200 Subject: [PATCH 01/12] [issue #45] Clean up imports statements. --- psydac/api/ast/evaluation.py | 84 ++++++++--------------------------- psydac/api/ast/expr.py | 69 ++++++++++------------------ psydac/api/ast/fem.py | 69 +++++++++++----------------- psydac/api/ast/glt.py | 66 ++++++++++----------------- psydac/api/ast/utilities.py | 39 ++++++++-------- psydac/api/printing/pycode.py | 9 ++-- 6 files changed, 113 insertions(+), 223 deletions(-) diff --git a/psydac/api/ast/evaluation.py b/psydac/api/ast/evaluation.py index e4f9aa931..f96a83c1b 100644 --- a/psydac/api/ast/evaluation.py +++ b/psydac/api/ast/evaluation.py @@ -1,78 +1,32 @@ -# TODO: - replace call to print_expression by SymbolicExpr (may need to use LogicalExpr) - from collections import OrderedDict -from itertools import groupby -import numpy as np - -from sympy import Basic -from sympy import symbols, Symbol, IndexedBase, Indexed, Function -from sympy import Mul, Add, Tuple, Min, Max, Pow -from sympy import Matrix, ImmutableDenseMatrix -from sympy import sqrt as sympy_sqrt -from sympy import S as sympy_S -from sympy import Integer, Float -from sympy.core.relational import Le, Ge -from sympy.logic.boolalg import And -from sympy import Mod, Abs -from sympy.core.function import AppliedUndef - -from pyccel.ast.core import Variable, IndexedVariable -from pyccel.ast.core import For -from pyccel.ast.core import Assign -from pyccel.ast.core import AugAssign -from pyccel.ast.core import Slice -from pyccel.ast.core import Range, Product -from pyccel.ast.core import FunctionDef -from pyccel.ast.core import FunctionCall -from pyccel.ast.core import Import -from pyccel.ast import Zeros -from pyccel.ast import Import -from pyccel.ast import DottedName -from pyccel.ast import Nil -from pyccel.ast import Len -from pyccel.ast import If, Is, Return -from pyccel.ast import String, Print, Shape -from pyccel.ast import Comment, NewLine -from pyccel.ast.core import _atomic -from pyccel.ast.utilities import build_types_decorator +from sympy import symbols +from sympy import Tuple, Pow -from sympde.core import Cross_3d -from sympde.core import Constant -from sympde.calculus import grad -from sympde.topology import Mapping -from sympde.topology import ScalarField -from sympde.topology import VectorField, IndexedVectorField -from sympde.topology import Boundary, BoundaryVector, NormalVector, TangentVector -from sympde.topology import Covariant, Contravariant -from sympde.topology import ElementArea -from sympde.topology import LogicalExpr -from sympde.topology import SymbolicExpr -from sympde.topology.derivatives import _partial_derivatives +from sympde.topology import Mapping +from sympde.topology import ScalarField +from sympde.topology import VectorField +from sympde.topology import SymbolicExpr +from sympde.topology.space import ScalarTestFunction from sympde.topology.derivatives import _logical_partial_derivatives -from sympde.topology.derivatives import get_max_partial_derivatives -from sympde.topology.space import ScalarFunctionSpace -from sympde.topology.space import ScalarTestFunction -from sympde.topology.space import VectorTestFunction -from sympde.topology.space import IndexedTestTrial -from sympde.topology.space import Trace -from sympde.topology.derivatives import print_expression -from sympde.topology.derivatives import get_atom_derivatives -from sympde.topology.derivatives import get_index_derivatives -from sympde.expr import BilinearForm, LinearForm, Functional, BasicForm - -from psydac.fem.splines import SplineSpace -from psydac.fem.tensor import TensorFemSpace -from psydac.fem.vector import ProductFemSpace - -from .basic import SplBasic +from sympde.topology.derivatives import _partial_derivatives # TODO: remove +from sympde.topology.derivatives import print_expression # TODO: remove + +from pyccel.ast.core import IndexedVariable +from pyccel.ast.core import For +from pyccel.ast.core import Assign +from pyccel.ast.core import Slice +from pyccel.ast.core import Range +from pyccel.ast.core import FunctionDef +from pyccel.ast.utilities import build_types_decorator + +from .basic import SplBasic from .utilities import build_pythran_types_header, variables from .utilities import filter_loops, filter_product from .utilities import rationalize_eval_mapping from .utilities import compute_atoms_expr_mapping from .utilities import compute_atoms_expr_field - #============================================================================== class EvalQuadratureMapping(SplBasic): diff --git a/psydac/api/ast/expr.py b/psydac/api/ast/expr.py index a8f2f3187..c23ac5f36 100644 --- a/psydac/api/ast/expr.py +++ b/psydac/api/ast/expr.py @@ -1,20 +1,9 @@ - -from collections import OrderedDict -from itertools import groupby -import string -import random -import numpy as np - -from sympy import Basic -from sympy import symbols, Symbol, IndexedBase, Indexed, Function -from sympy import Mul, Add, Tuple +from sympy import symbols, Symbol, IndexedBase +from sympy import Mul, Tuple from sympy import Matrix, ImmutableDenseMatrix -from sympy import sqrt as sympy_sqrt -from sympy import S as sympy_S -from sympy import simplify, expand from sympy.core.numbers import ImaginaryUnit -from pyccel.ast.core import Variable, IndexedVariable +from pyccel.ast.core import IndexedVariable from pyccel.ast.core import For from pyccel.ast.core import Assign from pyccel.ast.core import AugAssign @@ -22,42 +11,32 @@ from pyccel.ast.core import Range from pyccel.ast.core import FunctionDef from pyccel.ast.core import FunctionCall -from pyccel.ast import Zeros -from pyccel.ast import Import -from pyccel.ast import DottedName -from pyccel.ast import Nil -from pyccel.ast import Len -from pyccel.ast import If, Is, Return -from pyccel.ast import String, Print, Shape -from pyccel.ast import Comment, NewLine -from pyccel.ast.core import _atomic - -from sympde.core import Constant -from sympde.topology import ScalarField, VectorField -from sympde.topology import IndexedVectorField -from sympde.topology import Mapping -from sympde.expr import BilinearForm +from pyccel.ast import Import +from pyccel.ast import Nil +from pyccel.ast import Len +from pyccel.ast import If, Is, Return +from pyccel.ast.core import _atomic + +from sympde.core import Constant +from sympde.topology import ScalarField, VectorField +from sympde.topology import IndexedVectorField from sympde.topology.derivatives import _partial_derivatives from sympde.topology.derivatives import _logical_partial_derivatives from sympde.topology.derivatives import get_max_partial_derivatives -from sympde.topology.derivatives import print_expression from sympde.topology.derivatives import get_atom_derivatives from sympde.topology.derivatives import get_index_derivatives -from sympde.topology import LogicalExpr -from sympde.topology import SymbolicExpr -from sympde.topology import SymbolicDeterminant - -from .basic import SplBasic -from .utilities import random_string -from .utilities import build_pythran_types_header, variables -from .utilities import is_scalar_field, is_vector_field, is_mapping -from .utilities import math_atoms_as_str -#from .evaluation import EvalArrayVectorField -from .evaluation import EvalArrayMapping, EvalArrayField - -from psydac.fem.splines import SplineSpace -from psydac.fem.tensor import TensorFemSpace -from psydac.fem.vector import ProductFemSpace, VectorFemSpace +from sympde.topology.derivatives import print_expression # TODO: remove +from sympde.topology import LogicalExpr +from sympde.topology import SymbolicExpr +from sympde.topology import SymbolicDeterminant + +from .basic import SplBasic +from .utilities import random_string +from .utilities import build_pythran_types_header, variables +from .utilities import is_scalar_field, is_vector_field +from .utilities import math_atoms_as_str + +from psydac.fem.vector import ProductFemSpace #============================================================================== def compute_atoms_expr(atom, basis, indices, loc_indices, dim): diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index a3f7cd5ac..cec55180e 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -1,18 +1,10 @@ -# TODO: - replace call to print_expression by SymbolicExpr (may need to use LogicalExpr) - from collections import OrderedDict from itertools import groupby import numpy as np -from sympy import Basic -from sympy import symbols, Symbol, IndexedBase, Indexed, Function -from sympy import Mul, Add, Tuple, Min, Max, Pow +from sympy import symbols, Symbol, IndexedBase +from sympy import Mul, Tuple from sympy import Matrix, ImmutableDenseMatrix -from sympy import sqrt as sympy_sqrt -from sympy import S as sympy_S -from sympy import Integer, Float -from sympy.core.relational import Le, Ge -from sympy.logic.boolalg import And from sympy import Mod, Abs from sympy.core.function import AppliedUndef @@ -24,7 +16,6 @@ from pyccel.ast.core import Range, Product from pyccel.ast.core import FunctionDef from pyccel.ast.core import FunctionCall -from pyccel.ast.core import Import from pyccel.ast import Zeros from pyccel.ast import Import from pyccel.ast import DottedName @@ -32,49 +23,41 @@ from pyccel.ast import Len from pyccel.ast import If, Is, Return from pyccel.ast import String, Print, Shape -from pyccel.ast import Comment, NewLine +from pyccel.ast import Comment from pyccel.ast.core import _atomic from pyccel.ast.utilities import build_types_decorator -from sympde.core import Cross_3d -from sympde.core import Constant -from sympde.calculus import grad -from sympde.topology import Mapping -from sympde.topology import ScalarField -from sympde.topology import VectorField, IndexedVectorField -from sympde.topology import Boundary, BoundaryVector, NormalVector, TangentVector -from sympde.topology import Covariant, Contravariant -from sympde.topology import ElementArea -from sympde.topology import LogicalExpr -from sympde.topology import SymbolicExpr -from sympde.topology import UndefinedSpaceType +from sympde.core import Constant +from sympde.topology import ScalarField +from sympde.topology import VectorField, IndexedVectorField +from sympde.topology import Boundary, BoundaryVector, NormalVector, TangentVector +from sympde.topology import ElementArea +from sympde.topology import LogicalExpr +from sympde.topology import SymbolicExpr +from sympde.topology import UndefinedSpaceType +from sympde.topology.space import ScalarFunctionSpace, VectorFunctionSpace +from sympde.topology.space import ProductSpace +from sympde.topology.space import ScalarTestFunction +from sympde.topology.space import VectorTestFunction +from sympde.topology.space import IndexedTestTrial from sympde.topology.derivatives import _partial_derivatives -from sympde.topology.derivatives import _logical_partial_derivatives from sympde.topology.derivatives import get_max_partial_derivatives -from sympde.topology.space import ScalarFunctionSpace, VectorFunctionSpace -from sympde.topology.space import ProductSpace -from sympde.topology.space import ScalarTestFunction -from sympde.topology.space import VectorTestFunction -from sympde.topology.space import IndexedTestTrial -from sympde.topology.space import Trace -from sympde.topology.derivatives import print_expression -from sympde.topology.derivatives import get_atom_derivatives -from sympde.topology.derivatives import get_index_derivatives -from sympde.expr import BilinearForm, LinearForm, Functional, BasicForm +from sympde.topology.derivatives import print_expression # TODO: remove +from sympde.expr import BilinearForm, LinearForm, Functional from psydac.fem.splines import SplineSpace from psydac.fem.tensor import TensorFemSpace from psydac.fem.vector import ProductFemSpace -from .basic import SplBasic +from .basic import SplBasic from .evaluation import EvalQuadratureMapping, EvalQuadratureField, EvalQuadratureVectorField -from .utilities import random_string -from .utilities import build_pythran_types_header, variables -from .utilities import compute_normal_vector, compute_tangent_vector -from .utilities import select_loops, filter_product -from .utilities import compute_atoms_expr -from .utilities import is_scalar_field, is_vector_field -from .utilities import math_atoms_as_str +from .utilities import random_string +from .utilities import build_pythran_types_header, variables +from .utilities import compute_normal_vector, compute_tangent_vector +from .utilities import select_loops, filter_product +from .utilities import compute_atoms_expr +from .utilities import is_scalar_field, is_vector_field +from .utilities import math_atoms_as_str FunctionalForms = (BilinearForm, LinearForm, Functional) diff --git a/psydac/api/ast/glt.py b/psydac/api/ast/glt.py index a445c2463..d0baa64a8 100644 --- a/psydac/api/ast/glt.py +++ b/psydac/api/ast/glt.py @@ -1,69 +1,49 @@ - from collections import OrderedDict from itertools import groupby -import string -import random -import numpy as np -from sympy import Basic -from sympy import symbols, Symbol, IndexedBase, Indexed, Function -from sympy import Mul, Add, Tuple +from sympy import symbols, Symbol, IndexedBase +from sympy import Tuple from sympy import Matrix, ImmutableDenseMatrix -from sympy import sqrt as sympy_sqrt -from sympy import S as sympy_S from sympy import simplify, expand from sympy.core.numbers import ImaginaryUnit -from pyccel.ast.core import Variable, IndexedVariable +from pyccel.ast.core import IndexedVariable from pyccel.ast.core import For from pyccel.ast.core import Assign -from pyccel.ast.core import AugAssign from pyccel.ast.core import Slice from pyccel.ast.core import Range from pyccel.ast.core import FunctionDef from pyccel.ast.core import FunctionCall -from pyccel.ast import Zeros -from pyccel.ast import Import -from pyccel.ast import DottedName -from pyccel.ast import Nil -from pyccel.ast import Len -from pyccel.ast import If, Is, Return -from pyccel.ast import String, Print, Shape -from pyccel.ast import Comment, NewLine -from pyccel.ast.core import _atomic - -from sympde.core import Constant -from sympde.topology import ScalarField, VectorField -from sympde.topology import IndexedVectorField -from sympde.topology import Mapping -from sympde.expr import BilinearForm +from pyccel.ast import Zeros +from pyccel.ast import Import +from pyccel.ast import DottedName +from pyccel.ast import Nil +from pyccel.ast import Len +from pyccel.ast import If, Is, Return +from pyccel.ast.core import _atomic + +from sympde.topology import ScalarField, VectorField +from sympde.topology import IndexedVectorField from sympde.topology.derivatives import _partial_derivatives from sympde.topology.derivatives import _logical_partial_derivatives from sympde.topology.derivatives import get_max_partial_derivatives -from sympde.topology.derivatives import print_expression -from sympde.topology.derivatives import get_atom_derivatives -from sympde.topology.derivatives import get_index_derivatives -from sympde.topology import LogicalExpr -from sympde.topology import SymbolicExpr -from sympde.topology import SymbolicDeterminant - -#from gelato.core import gelatize +from sympde.topology.derivatives import print_expression # TODO: remove +from sympde.topology import LogicalExpr +from sympde.topology import SymbolicExpr +from sympde.topology import SymbolicDeterminant -from gelato.glt import BasicGlt from gelato.expr import gelatize -from .basic import SplBasic -from .utilities import random_string -from .utilities import build_pythran_types_header, variables -from .utilities import is_scalar_field, is_vector_field, is_mapping -from .utilities import math_atoms_as_str -#from .evaluation import EvalArrayVectorField +from .basic import SplBasic +from .utilities import random_string +from .utilities import build_pythran_types_header, variables +from .utilities import is_scalar_field, is_vector_field, is_mapping +from .utilities import math_atoms_as_str from .evaluation import EvalArrayMapping, EvalArrayField -from psydac.fem.splines import SplineSpace -from psydac.fem.tensor import TensorFemSpace from psydac.fem.vector import ProductFemSpace +#============================================================================== class GltKernel(SplBasic): def __new__(cls, expr, spaces, name=None, mapping=None, is_rational_mapping=None, backend=None): diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 4f0c2f3c5..283d34473 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -4,38 +4,35 @@ import numpy as np from sympy import Symbol, IndexedBase, Indexed -from sympy import Mul, Add, Tuple, Min, Max, Pow +from sympy import Mul, Tuple, Pow from sympy import Matrix from sympy import sqrt as sympy_sqrt from sympy.utilities.iterables import cartes +from sympde.topology.space import ScalarTestFunction +from sympde.topology.space import VectorTestFunction +from sympde.topology.space import IndexedTestTrial +from sympde.topology import ScalarField +from sympde.topology import VectorField, IndexedVectorField +from sympde.topology import Mapping +from sympde.topology.derivatives import _partial_derivatives +from sympde.topology.derivatives import _logical_partial_derivatives +from sympde.topology.derivatives import get_atom_derivatives +from sympde.topology.derivatives import get_index_derivatives +from sympde.topology.derivatives import print_expression # TODO: remove +from sympde.topology import LogicalExpr +from sympde.topology import SymbolicExpr +from sympde.core import Cross_3d + from pyccel.ast.core import Variable, IndexedVariable from pyccel.ast.core import For from pyccel.ast.core import Assign from pyccel.ast.core import AugAssign from pyccel.ast.core import Range, Product from pyccel.ast.core import _atomic -from pyccel.ast import Comment, NewLine - -from psydac.api.printing.pycode import pycode - - -from sympde.topology.derivatives import _partial_derivatives -from sympde.topology.derivatives import _logical_partial_derivatives -from sympde.topology.space import ScalarTestFunction -from sympde.topology.space import VectorTestFunction -from sympde.topology.space import IndexedTestTrial -from sympde.topology import ScalarField -from sympde.topology import VectorField, IndexedVectorField -from sympde.topology import Mapping -from sympde.topology.derivatives import print_expression -from sympde.topology.derivatives import get_atom_derivatives -from sympde.topology.derivatives import get_index_derivatives -from sympde.topology import LogicalExpr -from sympde.topology import SymbolicExpr -from sympde.core import Cross_3d - +from pyccel.ast import Comment +from psydac.api.printing.pycode import pycode #============================================================================== def random_string( n ): diff --git a/psydac/api/printing/pycode.py b/psydac/api/printing/pycode.py index fc102b80c..bcb8caf49 100644 --- a/psydac/api/printing/pycode.py +++ b/psydac/api/printing/pycode.py @@ -1,14 +1,11 @@ from sympy.core import Symbol -from sympy import Tuple from pyccel.codegen.printing.pycode import PythonCodePrinter as PyccelPythonCodePrinter -from sympde.calculus import Dot, Inner, Cross -from sympde.calculus import Grad, Rot, Curl, Div -from sympde.topology import Line, Square, Cube from sympde.topology.derivatives import _partial_derivatives -from sympde.topology.derivatives import print_expression +from sympde.topology.derivatives import print_expression # TODO: remove +#============================================================================== class PythonCodePrinter(PyccelPythonCodePrinter): def __init__(self, settings=None): @@ -117,7 +114,7 @@ def _print_AppliedUndef(self, expr): return '{fname}({args})'.format(fname=fname, args=args) # ......................................................... - +#============================================================================== def pycode(expr, **settings): """ Converts an expr to a string of Python code Parameters From 570737c4647b215ab20c967213f8809f917ee5ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 26 Aug 2019 18:46:32 +0200 Subject: [PATCH 02/12] [issue #45] Remove print_expression from pycode. --- psydac/api/printing/pycode.py | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/psydac/api/printing/pycode.py b/psydac/api/printing/pycode.py index bcb8caf49..685772010 100644 --- a/psydac/api/printing/pycode.py +++ b/psydac/api/printing/pycode.py @@ -3,7 +3,7 @@ from pyccel.codegen.printing.pycode import PythonCodePrinter as PyccelPythonCodePrinter from sympde.topology.derivatives import _partial_derivatives -from sympde.topology.derivatives import print_expression # TODO: remove +from sympde.topology import SymbolicExpr #============================================================================== class PythonCodePrinter(PyccelPythonCodePrinter): @@ -64,31 +64,25 @@ def _print_GltInterface(self, expr): def _print_dx(self, expr): arg = expr.args[0] if isinstance(arg, _partial_derivatives): - arg = print_expression(arg, mapping_name=False) - + arg = SymbolicExpr(arg).name else: arg = self._print(arg) + '_' - return arg + 'x' def _print_dy(self, expr): arg = expr.args[0] if isinstance(arg, _partial_derivatives): - arg = print_expression(arg, mapping_name=False) - + arg = SymbolicExpr(arg).name else: arg = self._print(arg) + '_' - return arg + 'y' def _print_dz(self, expr): arg = expr.args[0] if isinstance(arg, _partial_derivatives): - arg = print_expression(arg, mapping_name=False) - + arg = SymbolicExpr(arg).name else: arg = self._print(arg) + '_' - return arg + 'z' def _print_IndexedTestTrial(self, expr): From 0f46ed877324b12f3830ee4ee70a1b5a1d85c0e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 28 Aug 2019 15:08:14 +0200 Subject: [PATCH 03/12] [issue #45] Remove print_expression from 'api.ast.utilities'. --- psydac/api/ast/utilities.py | 132 +++++++++++++++++++----------------- 1 file changed, 69 insertions(+), 63 deletions(-) diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 283d34473..56765ff3d 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -19,7 +19,8 @@ from sympde.topology.derivatives import _logical_partial_derivatives from sympde.topology.derivatives import get_atom_derivatives from sympde.topology.derivatives import get_index_derivatives -from sympde.topology.derivatives import print_expression # TODO: remove +from sympde.topology.derivatives import get_atom_logical_derivatives +from sympde.topology.derivatives import get_index_logical_derivatives from sympde.topology import LogicalExpr from sympde.topology import SymbolicExpr from sympde.core import Cross_3d @@ -32,7 +33,7 @@ from pyccel.ast.core import _atomic from pyccel.ast import Comment -from psydac.api.printing.pycode import pycode +from psydac.api.printing.pycode import pycode # TODO: remove #============================================================================== def random_string( n ): @@ -176,20 +177,23 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test, map_stmts = [] if mapping: for atom in atomic_exprs: + if isinstance(atom, _partial_derivatives): - name = print_expression(atom) - rhs = LogicalExpr(mapping, atom) + lhs = SymbolicExpr(atom) + rhs_p = LogicalExpr(mapping, atom) + # we look for new_atoms that must be added to atomic_exprs # because we need them in the maps stmts - logical_atoms = _atomic(rhs, cls=_logical_partial_derivatives) - rhs = SymbolicExpr(rhs) - map_stmts += [Assign(Symbol(name), rhs)] - for atom in logical_atoms: - ls = _atomic(atom, Symbol) + logical_atoms = _atomic(rhs_p, cls=_logical_partial_derivatives) + for a in logical_atoms: + ls = _atomic(a, Symbol) assert len(ls) == 1 if isinstance(ls[0], cls): - new_atoms += [logical2physical(atom)] + new_atoms += [logical2physical(a)] + rhs = SymbolicExpr(rhs_p) + map_stmts += [Assign(lhs, rhs)] + atomic_exprs = {*atomic_exprs, *new_atoms} assigns = [] @@ -209,9 +213,13 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test, args = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)] # ... assign basis on quad point - logical = not( mapping is None ) - name = print_expression(atom, logical=logical) - assigns += [Assign(Symbol(name), Mul(*args))] + if mapping: + subs = dict(zip(_partial_derivatives, _logical_partial_derivatives)) + atom = atom.subs(subs) + + lhs = SymbolicExpr(atom) + rhs = Mul(*args) + assigns += [Assign(lhs, rhs)] # ... return assigns, map_stmts @@ -275,51 +283,52 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad, # ... map basis function new_atoms = [] for atom in atomic_exprs: + if isinstance(atom, _partial_derivatives): - name = print_expression(atom) - rhs = LogicalExpr(mapping, atom) + lhs = SymbolicExpr(atom) + rhs_p = LogicalExpr(mapping, atom) - logical_atoms = _atomic(rhs, cls=_logical_partial_derivatives) - rhs = SymbolicExpr(rhs) - sym = Symbol(name) - map_stmts += [Assign(sym, rhs)] + logical_atoms = _atomic(rhs_p, cls=_logical_partial_derivatives) + for a in logical_atoms: + ls = _atomic(a, Symbol) + if isinstance(ls[0], cls): + new_atoms += [logical2physical(a)] - if not pycode(atom) == name: - var = Symbol(pycode(atom)) - map_stmts += [Assign(var, sym)] + rhs = SymbolicExpr(rhs_p) + map_stmts += [Assign(lhs, rhs)] - for atom in logical_atoms: - ls = _atomic(atom, Symbol) - if isinstance(ls[0], cls): - new_atoms += [logical2physical(atom)] + # Temporary hack: add 'phi_yx = phi_xy' where necessary + # TODO: remove pycode calls from Psydac's AST + if pycode(atom) != lhs.name: + var = Symbol(pycode(atom)) + map_stmts += [Assign(var, lhs)] atomic_exprs = {*atomic_exprs, *new_atoms} - atomic_exprs = {print_expression(a):a for a in atomic_exprs} + # Make sure that we only pick one between 'dx(dy(u))' and 'dy(dx(u))' + atomic_exprs = {SymbolicExpr(a).name : a for a in atomic_exprs} atomic_exprs = tuple(atomic_exprs.values()) for atom in atomic_exprs: if is_scalar_field(atom): - field = list(atom.atoms(ScalarField))[0] - field_name = 'coeff_' + print_expression(field) + field = atom.atoms(ScalarField).pop() + field_name = 'coeff_' + SymbolicExpr(field).name base = field elif is_vector_field(atom): - field = list(atom.atoms(IndexedVectorField))[0] - field_name = 'coeff_' + print_expression(field) + field = atom.atoms(IndexedVectorField).pop() + field_name = 'coeff_' + SymbolicExpr(field).name base = field.base else: raise TypeError('atom must be either scalar or vector field') # ... - - test_fun = atom.subs(base, test_function) - name = print_expression(test_fun) - test_fun = Symbol(name) + test_fun = SymbolicExpr(atom.subs(base, test_function)) # ... # ... + # TODO: verify if 'get_index_logical_derivatives' should be used instead orders = [*get_index_derivatives(atom).values()] args = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)] inits += [Assign(test_fun, Mul(*args))] @@ -327,7 +336,7 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad, # ... args = [IndexedBase(field_name)[idxs], test_fun] - val_name = print_expression(atom) + '_values' + val_name = SymbolicExpr(atom).name + '_values' val = IndexedBase(val_name)[indices_quad] updates += [AugAssign(val,'+',Mul(*args))] # ... @@ -371,30 +380,29 @@ def compute_atoms_expr_mapping(atomic_exprs, indices_quad, """ - _print = lambda i: print_expression(i, mapping_name=False) inits = [] updates = [] for atom in atomic_exprs: - element = get_atom_derivatives(atom) - element_name = 'coeff_' + _print(element) + + element = get_atom_logical_derivatives(atom) + element_name = 'coeff_' + SymbolicExpr(element).name # ... test_fun = atom.subs(element, test_function) - name = print_expression(test_fun, logical=True) - test_fun = Symbol(name) + test_fun = SymbolicExpr(test_fun) # ... # ... - orders = [*get_index_derivatives(atom).values()] + orders = [*get_index_logical_derivatives(atom).values()] args = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)] inits += [Assign(test_fun, Mul(*args))] # ... # ... - args = [IndexedBase(element_name)[idxs], test_fun] - val_name = _print(atom) + '_values' + val_name = SymbolicExpr(atom).name + '_values' val = IndexedBase(val_name)[indices_quad] - updates += [AugAssign(val,'+',Mul(*args))] + expr = IndexedBase(element_name)[idxs] * test_fun + updates += [AugAssign(val, '+', expr)] # ... return inits, updates @@ -402,11 +410,9 @@ def compute_atoms_expr_mapping(atomic_exprs, indices_quad, #============================================================================== def rationalize_eval_mapping(mapping, nderiv, space, indices_quad): - _print = lambda i: print_expression(i, mapping_name=False) - M = mapping dim = space.ldim - ops = _partial_derivatives[:dim] + ops = _logical_partial_derivatives[:dim] # ... mapping components and their derivatives components = [M[i] for i in range(0, dim)] @@ -444,7 +450,7 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad): # declarations stmts += [Comment('declarations')] for atom in elements + weights_elements: - atom_name = _print(atom) + atom_name = SymbolicExpr(atom).name val_name = atom_name + '_values' val = IndexedBase(val_name)[indices_quad] @@ -455,9 +461,9 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad): stmts += [Comment('rationalize')] # 0 order terms - w = Symbol( _print( weights ) ) for i in range(dim): - u = Symbol( _print( M[i] ) ) + w = SymbolicExpr(weights) + u = SymbolicExpr(M[i]) val_name = u.name + '_values' val = IndexedBase(val_name)[indices_quad] @@ -468,12 +474,12 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad): # 1 order terms if nderiv >= 1: for d in ops: - w = Symbol( _print( weights ) ) - dw = Symbol( _print( d(weights) ) ) + w = SymbolicExpr( weights ) + dw = SymbolicExpr(d(weights)) for i in range(dim): - u = Symbol( _print( M[i] ) ) - du = Symbol( _print( d(M[i]) ) ) + u = SymbolicExpr( M[i] ) + du = SymbolicExpr(d(M[i])) val_name = du.name + '_values' val = IndexedBase(val_name)[indices_quad] @@ -485,16 +491,16 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad): if nderiv >= 2: for e, d1 in enumerate(ops): for d2 in ops[:e+1]: - w = Symbol( _print( weights ) ) - d1w = Symbol( _print( d1(weights) ) ) - d2w = Symbol( _print( d2(weights) ) ) - d1d2w = Symbol( _print( d1(d2(weights)) ) ) + w = SymbolicExpr( weights ) + d1w = SymbolicExpr( d1(weights) ) + d2w = SymbolicExpr( d2(weights) ) + d1d2w = SymbolicExpr(d1(d2(weights))) for i in range(dim): - u = Symbol( _print( M[i] ) ) - d1u = Symbol( _print( d1(M[i]) ) ) - d2u = Symbol( _print( d2(M[i]) ) ) - d1d2u = Symbol( _print( d1(d2(M[i])) ) ) + u = SymbolicExpr( M[i] ) + d1u = SymbolicExpr( d1(M[i]) ) + d2u = SymbolicExpr( d2(M[i]) ) + d1d2u = SymbolicExpr(d1(d2(M[i]))) val_name = d1d2u.name + '_values' val = IndexedBase(val_name)[indices_quad] From 4b592041c9f9bb27e4610f5099b194b9ccd65b56 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 28 Aug 2019 15:44:42 +0200 Subject: [PATCH 04/12] [issue #45] Remove print_expression from 'api.ast.evaluation'. --- psydac/api/ast/evaluation.py | 51 ++++++++++++++++-------------------- 1 file changed, 23 insertions(+), 28 deletions(-) diff --git a/psydac/api/ast/evaluation.py b/psydac/api/ast/evaluation.py index f96a83c1b..e70012448 100644 --- a/psydac/api/ast/evaluation.py +++ b/psydac/api/ast/evaluation.py @@ -9,8 +9,6 @@ from sympde.topology import SymbolicExpr from sympde.topology.space import ScalarTestFunction from sympde.topology.derivatives import _logical_partial_derivatives -from sympde.topology.derivatives import _partial_derivatives # TODO: remove -from sympde.topology.derivatives import print_expression # TODO: remove from pyccel.ast.core import IndexedVariable from pyccel.ast.core import For @@ -54,7 +52,7 @@ def __new__(cls, space, mapping, discrete_boundary=None, name=None, # ... # ... - ops = _partial_derivatives[:dim] + ops = _logical_partial_derivatives[:dim] M = mapping components = [M[i] for i in range(0, dim)] @@ -150,10 +148,8 @@ def _initialize(self): space = self.space dim = space.ldim - _print = lambda i: print_expression(i, mapping_name=False) - mapping_atoms = [_print(f) for f in self.components] - mapping_str = [_print(f) for f in self.elements] -# mapping_str = sorted([_print(f) for f in self.elements]) + mapping_atoms = [SymbolicExpr(f).name for f in self.components] + mapping_str = [SymbolicExpr(f).name for f in self.elements ] # ... declarations degrees = variables( 'p1:%s'%(dim+1), 'int') @@ -183,7 +179,7 @@ def _initialize(self): # ... nderiv = self.nderiv - ops = _partial_derivatives[:dim] + ops = _logical_partial_derivatives[:dim] if nderiv > 0: weights_elements += [d(weights_pts) for d in ops] @@ -196,7 +192,7 @@ def _initialize(self): raise NotImplementedError('TODO') # ... - mapping_weights_str = [_print(f) for f in weights_elements] + mapping_weights_str = [SymbolicExpr(f).name for f in weights_elements] mapping_wvalues = variables(['{}_values'.format(f) for f in mapping_weights_str], dtype='real', rank=dim, cls=IndexedVariable) @@ -221,11 +217,12 @@ def _initialize(self): # ... # ... - Nj = ScalarTestFunction(space, name='Nj') - body = [] - init_basis = OrderedDict() - updates = [] + Nj = ScalarTestFunction(space, name='Nj') + body = [] + init_basis = OrderedDict() + updates = [] atomic_exprs = self.elements + weights_elements + inits, updates = compute_atoms_expr_mapping(atomic_exprs, indices_quad, indices_basis, basis, Nj) @@ -310,7 +307,6 @@ def __new__(cls, space, fields, discrete_boundary=None, name=None, obj._backend = backend obj._func = obj._initialize() - return obj @property @@ -385,8 +381,8 @@ def _initialize(self): self._fields = Tuple(*fields) - fields_str = [print_expression(f) for f in self._fields] - fields_val = variables(['{}_values'.format(f) for f in fields_str], + fields_str = [SymbolicExpr(f).name for f in self._fields] + fields_val = variables(['{}_values'.format(f) for f in fields_str], dtype='real', rank=dim, cls=IndexedVariable) @@ -493,7 +489,7 @@ def _initialize(self): indices_quad = variables('g1:%s'%(dim+1), 'int') basis = variables('basis1:%s'%(dim+1), dtype='real', rank=3, cls=IndexedVariable) - coeffs = ['coeff_{}'.format(print_expression(f)) for f in vector_field_atoms] + coeffs = ['coeff_{}'.format(SymbolicExpr(f).name) for f in vector_field_atoms] vector_fields_coeffs = variables(coeffs, dtype='real', rank=dim, cls=IndexedVariable) # ... @@ -522,7 +518,7 @@ def _initialize(self): init_map[str(stmt.lhs)] = stmt self._vector_fields = Tuple(*vector_fields) - vector_fields_str = [print_expression(f) for f in self.vector_fields] + vector_fields_str = [SymbolicExpr(f).name for f in self.vector_fields] vector_fields_val = variables(['{}_values'.format(f) for f in vector_fields_str], dtype='real', rank=dim, cls=IndexedVariable) @@ -580,6 +576,7 @@ def _create_loop(indices, ranges, body): return body #============================================================================== +# NOTE: this is used in module 'psydac.api.ast.glt' class EvalArrayField(SplBasic): def __new__(cls, space, fields, discrete_boundary=None, name=None, @@ -599,7 +596,6 @@ def __new__(cls, space, fields, discrete_boundary=None, name=None, obj._backend = backend obj._func = obj._initialize() - return obj @property @@ -634,7 +630,7 @@ def _initialize(self): mapping = self.mapping field_atoms = self.fields.atoms(ScalarField) - fields_str = sorted([print_expression(f) for f in self.fields]) + fields_str = sorted([SymbolicExpr(f).name for f in self.fields]) # ... declarations degrees = variables( 'p1:%s'%(dim+1), 'int') @@ -725,6 +721,7 @@ def _initialize(self): #============================================================================== +# NOTE: this is used in module 'psydac.api.ast.glt' class EvalArrayMapping(SplBasic): def __new__(cls, space, mapping, name=None, @@ -749,7 +746,7 @@ def __new__(cls, space, mapping, name=None, # ... # ... - ops = _partial_derivatives[:dim] + ops = _logical_partial_derivatives[:dim] M = mapping components = [M[i] for i in range(0, dim)] @@ -834,10 +831,8 @@ def _initialize(self): space = self.space dim = space.ldim - _print = lambda i: print_expression(i, mapping_name=False) - mapping_atoms = [_print(f) for f in self.components] - mapping_str = [_print(f) for f in self.elements] -# mapping_str = sorted([_print(f) for f in self.elements]) + mapping_atoms = [SymbolicExpr(f).name for f in self.components] + mapping_str = [SymbolicExpr(f).name for f in self.elements ] # ... declarations degrees = variables( 'p1:%s'%(dim+1), 'int') @@ -871,7 +866,7 @@ def _initialize(self): # ... nderiv = self.nderiv - ops = _partial_derivatives[:dim] + ops = _logical_partial_derivatives[:dim] if nderiv > 0: weights_elements += [d(weights_pts) for d in ops] @@ -884,7 +879,7 @@ def _initialize(self): raise NotImplementedError('TODO') # ... - mapping_weights_str = [_print(f) for f in weights_elements] + mapping_weights_str = [SymbolicExpr(f).name for f in weights_elements] mapping_wvalues = variables(['{}_values'.format(f) for f in mapping_weights_str], dtype='real', rank=dim, cls=IndexedVariable) @@ -937,8 +932,8 @@ def _initialize(self): assign_spans = [] for x, i_span, span in zip(indices_quad, i_spans, spans): assign_spans += [Assign(i_span, span[x])] - body = assign_spans + body + # put the body in for loops of quadrature points body = _create_loop(indices_quad, ranges_quad, body) From 8605869d8166e71974a294c86e33c1915e3f5796 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 28 Aug 2019 15:55:54 +0200 Subject: [PATCH 05/12] [issue #45] Update kernel AST when mapping is present. --- psydac/api/ast/fem.py | 21 +++++++-------------- psydac/api/ast/glt.py | 21 +++++++-------------- 2 files changed, 14 insertions(+), 28 deletions(-) diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index cec55180e..2d6b5e790 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -751,21 +751,14 @@ def _initialize(self): # ... # ... - mapping_elements = () - mapping_coeffs = () - mapping_values = () if mapping: - _eval = self.eval_mapping - _print = lambda i: print_expression(i, mapping_name=False) - - mapping_elements = [_print(i) for i in _eval.elements] - mapping_elements = symbols(tuple(mapping_elements)) - - mapping_coeffs = [_print(i) for i in _eval.mapping_coeffs] - mapping_coeffs = variables(mapping_coeffs, dtype='real', rank=dim, cls=IndexedVariable) - - mapping_values = [_print(i) for i in _eval.mapping_values] - mapping_values = variables(mapping_values, dtype='real', rank=dim, cls=IndexedVariable) + mapping_elements = [SymbolicExpr(i) for i in self.eval_mapping.elements] + mapping_coeffs = self.eval_mapping.mapping_coeffs + mapping_values = self.eval_mapping.mapping_values + else: + mapping_elements = () + mapping_coeffs = () + mapping_values = () # ... # ... diff --git a/psydac/api/ast/glt.py b/psydac/api/ast/glt.py index d0baa64a8..d9fc17c13 100644 --- a/psydac/api/ast/glt.py +++ b/psydac/api/ast/glt.py @@ -369,21 +369,14 @@ def _initialize(self): # ... # ... - mapping_elements = () - mapping_coeffs = () - mapping_values = () if mapping: - _eval = self.eval_mapping - _print = lambda i: print_expression(i, mapping_name=False) - - mapping_elements = [_print(i) for i in _eval.elements] - mapping_elements = symbols(tuple(mapping_elements)) - - mapping_coeffs = [_print(i) for i in _eval.mapping_coeffs] - mapping_coeffs = variables(mapping_coeffs, dtype='real', rank=dim, cls=IndexedVariable) - - mapping_values = [_print(i) for i in _eval.mapping_values] - mapping_values = variables(mapping_values, dtype='real', rank=dim, cls=IndexedVariable) + mapping_elements = [SymbolicExpr(i) for i in self.eval_mapping.elements] + mapping_coeffs = self.eval_mapping.mapping_coeffs + mapping_values = self.eval_mapping.mapping_values + else: + mapping_elements = () + mapping_coeffs = () + mapping_values = () # ... # # ... From 0457a0711a2d6967bf3b7549b5146e2f092c382c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 28 Aug 2019 18:43:44 +0200 Subject: [PATCH 06/12] [issue #45] Cleanup and bugfix in 'compute_atoms_expr'. --- psydac/api/ast/utilities.py | 65 +++++++++++++++++++++---------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 56765ff3d..03d9ceef5 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -148,8 +148,8 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test, Returns ------- - assignments : - list of assignments of the atomic expression evaluated in the quadature points + inits : + list of assignments of the atomic expression evaluated in the quadrature points map_stmts : list of assigments of atomic expression in case of mapping @@ -169,16 +169,23 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test, for atom in atomic_exprs: if not isinstance(atom, cls): raise TypeError('atom must be of type {}'.format(str(cls))) - # ... - atomic_exprs = list(atomic_exprs) - new_atoms = [] - # ... map basis function - map_stmts = [] + # If there is a mapping, compute [dx(u), dy(u), dz(u)] as functions + # of [dx1(u), dx2(u), dx3(u)], and store results into intermediate + # variables [u_x, u_y, u_z]. (Same thing is done for higher derivatives.) + # + # Accordingly, we create a new list of atoms where all partial derivatives + # are taken with respect to the logical coordinates. if mapping: + + new_atoms = set() + map_stmts = [] + get_index = get_index_logical_derivatives + get_atom = get_atom_logical_derivatives + for atom in atomic_exprs: - if isinstance(atom, _partial_derivatives): + if isinstance(atom, _partial_derivatives): lhs = SymbolicExpr(atom) rhs_p = LogicalExpr(mapping, atom) @@ -189,18 +196,26 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test, ls = _atomic(a, Symbol) assert len(ls) == 1 if isinstance(ls[0], cls): - new_atoms += [logical2physical(a)] + new_atoms.add(a) rhs = SymbolicExpr(rhs_p) map_stmts += [Assign(lhs, rhs)] - atomic_exprs = {*atomic_exprs, *new_atoms} - - assigns = [] - for atom in atomic_exprs: + else: + new_atoms.add(atom) - orders = [*get_index_derivatives(atom).values()] - a = get_atom_derivatives(atom) + else: + new_atoms = atomic_exprs + map_stmts = [] + get_index = get_index_derivatives + get_atom = get_atom_derivatives + + # Create a list of statements for initialization of the point values, + # for each of the atoms in our (possibly new) list. + inits = [] + for atom in new_atoms: + orders = [*get_index(atom).values()] + a = get_atom(atom) test = _get_name(a) in [_get_name(f) for f in test_function] if test or is_linear: @@ -210,20 +225,14 @@ def compute_atoms_expr(atomic_exprs, indices_quad, indices_test, basis = basis_trial idxs = indices_trial - args = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)] - - # ... assign basis on quad point - if mapping: - subs = dict(zip(_partial_derivatives, _logical_partial_derivatives)) - atom = atom.subs(subs) - - lhs = SymbolicExpr(atom) - rhs = Mul(*args) - assigns += [Assign(lhs, rhs)] - - # ... - return assigns, map_stmts + args = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)] + lhs = SymbolicExpr(atom) + rhs = Mul(*args) + inits += [Assign(lhs, rhs)] + # Return the initialization statements, and the additional initialization + # of intermediate variables in case of mapping + return inits, map_stmts #============================================================================== def compute_atoms_expr_field(atomic_exprs, indices_quad, From 22f374f764176d9c591503c69c38b5c98434786e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 28 Aug 2019 18:58:50 +0200 Subject: [PATCH 07/12] [issue #45] Remove print_expression from 'api.ast.fem'. --- psydac/api/ast/fem.py | 65 +++++++++++++++---------------------------- 1 file changed, 23 insertions(+), 42 deletions(-) diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index 2d6b5e790..f7443b3e6 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -41,8 +41,8 @@ from sympde.topology.space import VectorTestFunction from sympde.topology.space import IndexedTestTrial from sympde.topology.derivatives import _partial_derivatives +from sympde.topology.derivatives import _logical_partial_derivatives from sympde.topology.derivatives import get_max_partial_derivatives -from sympde.topology.derivatives import print_expression # TODO: remove from sympde.expr import BilinearForm, LinearForm, Functional from psydac.fem.splines import SplineSpace @@ -181,45 +181,26 @@ def init_loop_support(indices_elm, n_elements, # TODO take exponent to 1/dim def area_eval_mapping(mapping, area, dim, indices_quad, weight): - _print = lambda i: print_expression(i, mapping_name=False) - - M = mapping - ops = _partial_derivatives[:dim] + stmts = [] - # ... mapping components and their derivatives - elements = [d(M[i]) for d in ops for i in range(0, dim)] - # ... + # mapping components and their derivatives + ops = _logical_partial_derivatives[:dim] + elements = [d(mapping[i]) for d in ops for i in range(0, dim)] - stmts = [] # declarations stmts += [Comment('declarations')] - for atom in elements: - atom_name = _print(atom) - val_name = atom_name + '_values' - val = IndexedBase(val_name)[indices_quad] - - stmt = Assign(atom_name, val) - stmts += [stmt] - - # ... inv jacobian - jac = mapping.det_jacobian - rdim = mapping.rdim - ops = _partial_derivatives[:rdim] - elements = [d(mapping[i]) for d in ops for i in range(0, rdim)] for e in elements: - new = print_expression(e, mapping_name=False) - new = Symbol(new) - jac = jac.subs(e, new) - # ... + lhs = SymbolicExpr(e) + rhs_name = lhs.name + '_values' + rhs = IndexedBase(rhs_name)[indices_quad] + stmts += [Assign(lhs, rhs)] - # ... + # jacobian determinant + jac = SymbolicExpr(mapping.det_jacobian) stmts += [AugAssign(area, '+', Abs(jac) * weight)] - # ... return stmts - - #============================================================================== # target is used when there are multiple expression (domain/boundaries) class Kernel(SplBasic): @@ -550,21 +531,19 @@ def _initialize(self): # update dependencies self._dependencies += self.eval_fields - # TODO use print_expression - - - fields_str = tuple(map(print_expression, fields)) - fields_logical_str = [print_expression(f, logical=True) for f in - fields] - # ... + # TODO: remove these? + d_subs = dict(zip(_partial_derivatives, _logical_partial_derivatives)) + fields_logical = tuple(f.subs(d_subs) for f in fields) + fields_str = tuple(SymbolicExpr(f).name for f in fields) + fields_logical_str = tuple(SymbolicExpr(f).name for f in fields_logical) # ... vector_field_atoms = tuple(expr.atoms(VectorField)) vector_fields = [] + # ... create EvalQuadratureVectorField self._eval_vector_fields = [] - if atomic_expr_vector_field: keyfunc = lambda F: F.space.name data = sorted(vector_field_atoms, key=keyfunc) @@ -589,9 +568,11 @@ def _initialize(self): # update dependencies self._dependencies += self.eval_vector_fields - vector_fields_str = tuple(print_expression(i) for i in vector_fields) - vector_fields_logical_str = [print_expression(f, logical=True) for f in - vector_fields] + + # TODO: remove these? + vector_fields_logical = tuple(f.subs(d_subs) for f in vector_fields) + vector_fields_str = tuple(SymbolicExpr(f).name for f in vector_fields) + vector_fields_logical_str = tuple(SymbolicExpr(f).name for f in vector_fields_logical) # ... # ... TODO add it as a method to basic class @@ -706,7 +687,7 @@ def _initialize(self): vector_fields_logical = symbols(vector_fields_logical_str) vector_field_atoms = [f[i] for f in vector_field_atoms for i in range(0, dim)] - coeffs = ['coeff_{}'.format(print_expression(f)) for f in vector_field_atoms] + coeffs = ['coeff_{}'.format(SymbolicExpr(f).name) for f in vector_field_atoms] vector_fields_coeffs = variables(coeffs, dtype='real', rank=dim, cls=IndexedVariable) vector_fields_val = variables(['{}_values'.format(f) for f in vector_fields_str], From e4678886eab406e84568b192fecdf1d052e88aa3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 28 Aug 2019 19:03:00 +0200 Subject: [PATCH 08/12] [issue #45] Remove print_expression from 'api.ast.expr'. --- psydac/api/ast/expr.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/psydac/api/ast/expr.py b/psydac/api/ast/expr.py index c23ac5f36..a8c913a04 100644 --- a/psydac/api/ast/expr.py +++ b/psydac/api/ast/expr.py @@ -25,7 +25,6 @@ from sympde.topology.derivatives import get_max_partial_derivatives from sympde.topology.derivatives import get_atom_derivatives from sympde.topology.derivatives import get_index_derivatives -from sympde.topology.derivatives import print_expression # TODO: remove from sympde.topology import LogicalExpr from sympde.topology import SymbolicExpr from sympde.topology import SymbolicDeterminant @@ -247,8 +246,8 @@ def _initialize(self): self._fields = tuple(expr.atoms(ScalarField)) self._vector_fields = tuple(expr.atoms(VectorField)) # ... - fields_str = tuple(map(print_expression, atomic_expr_field)) - vector_fields_str = tuple(map(print_expression, atomic_expr_vector_field)) + fields_str = tuple(SymbolicExpr(f).name for f in atomic_expr_field) + vector_fields_str = tuple(SymbolicExpr(f).name for f in atomic_expr_vector_field) fields = symbols(fields_str) vector_fields = symbols(vector_fields_str) From 6640df009db34e0f7aafabd6871fc06aaff8ec16 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 28 Aug 2019 19:14:18 +0200 Subject: [PATCH 09/12] [issue #45] Remove print_expression from 'api.ast.glt'. --- psydac/api/ast/glt.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/psydac/api/ast/glt.py b/psydac/api/ast/glt.py index d9fc17c13..f24251c5c 100644 --- a/psydac/api/ast/glt.py +++ b/psydac/api/ast/glt.py @@ -27,7 +27,6 @@ from sympde.topology.derivatives import _partial_derivatives from sympde.topology.derivatives import _logical_partial_derivatives from sympde.topology.derivatives import get_max_partial_derivatives -from sympde.topology.derivatives import print_expression # TODO: remove from sympde.topology import LogicalExpr from sympde.topology import SymbolicExpr from sympde.topology import SymbolicDeterminant @@ -280,10 +279,11 @@ def _initialize(self): # ... # ... - fields_str = sorted(tuple(map(print_expression, atomic_expr_field))) - fields_logical_str = sorted([print_expression(f, logical=True) for f in - atomic_expr_field]) - field_atoms = tuple(expr.atoms(ScalarField)) + d_subs = dict(zip(_partial_derivatives, _logical_partial_derivatives)) + atomic_expr_field_logical = tuple(f.subs(d_subs) for f in atomic_expr_field) + fields_str = tuple(sorted(SymbolicExpr(f).name for f in atomic_expr_field)) + fields_logical_str = tuple(sorted(SymbolicExpr(f).name for f in atomic_expr_field_logical)) + field_atoms = tuple(expr.atoms(ScalarField)) # ... # ... create EvalArrayField @@ -361,7 +361,7 @@ def _initialize(self): # vector_fields_logical = symbols(vector_fields_logical_str) # # vector_field_atoms = [f[i] for f in vector_field_atoms for i in range(0, dim)] -# coeffs = ['coeff_{}'.format(print_expression(f)) for f in vector_field_atoms] +# coeffs = ['coeff_{}'.format(SymbolicExpr(f).name) for f in vector_field_atoms] # vector_fields_coeffs = variables(coeffs, dtype='real', rank=dim, cls=IndexedVariable) # # vector_fields_val = variables(['{}_values'.format(f) for f in vector_fields_str], From be5903f02da35e7fedb228db42fed58abd85a1f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 29 Aug 2019 19:05:47 +0200 Subject: [PATCH 10/12] [issue #45] Cleanup and bugfix in 'compute_atoms_expr_field'. --- psydac/api/ast/utilities.py | 69 ++++++++++++++++++++++++------------- 1 file changed, 45 insertions(+), 24 deletions(-) diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 03d9ceef5..b0575685d 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -62,6 +62,9 @@ def is_scalar_field(expr): if isinstance(expr, _partial_derivatives): return is_scalar_field(expr.args[0]) + elif isinstance(expr, _logical_partial_derivatives): + return is_scalar_field(expr.args[0]) + elif isinstance(expr, ScalarField): return True @@ -73,6 +76,9 @@ def is_vector_field(expr): if isinstance(expr, _partial_derivatives): return is_vector_field(expr.args[0]) + elif isinstance(expr, _logical_partial_derivatives): + return is_vector_field(expr.args[0]) + elif isinstance(expr, (VectorField, IndexedVectorField)): return True @@ -273,10 +279,10 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad, list of augmented assignments which are updated in each loop iteration map_stmts : - list of assigments of atomic expression in case of mapping + list of assignments of atomic expression in case of mapping - atomic_exprs: - list of new atomic expressions that were introduced in case of a mapping + new_atoms: + updated list of atomic expressions (some were introduced in case of a mapping) """ inits = [] @@ -288,57 +294,72 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad, IndexedVectorField, VectorField) + # If there is a mapping, compute [dx(u), dy(u), dz(u)] as functions + # of [dx1(u), dx2(u), dx3(u)], and store results into intermediate + # variables [u_x, u_y, u_z]. (Same thing is done for higher derivatives.) + # + # Accordingly, we create a new list of atoms where all partial derivatives + # are taken with respect to the logical coordinates. if mapping: - # ... map basis function - new_atoms = [] + + new_atoms = set() + map_stmts = [] + get_index = get_index_logical_derivatives + get_atom = get_atom_logical_derivatives + for atom in atomic_exprs: if isinstance(atom, _partial_derivatives): lhs = SymbolicExpr(atom) rhs_p = LogicalExpr(mapping, atom) + # we look for new_atoms that must be added to atomic_exprs + # because we need them in the maps stmts logical_atoms = _atomic(rhs_p, cls=_logical_partial_derivatives) for a in logical_atoms: ls = _atomic(a, Symbol) + assert len(ls) == 1 if isinstance(ls[0], cls): - new_atoms += [logical2physical(a)] + new_atoms.add(a) rhs = SymbolicExpr(rhs_p) map_stmts += [Assign(lhs, rhs)] - # Temporary hack: add 'phi_yx = phi_xy' where necessary - # TODO: remove pycode calls from Psydac's AST - if pycode(atom) != lhs.name: - var = Symbol(pycode(atom)) - map_stmts += [Assign(var, lhs)] + else: + new_atoms.add(atom) - atomic_exprs = {*atomic_exprs, *new_atoms} + else: + new_atoms = atomic_exprs + map_stmts = [] + get_index = get_index_derivatives + get_atom = get_atom_derivatives - # Make sure that we only pick one between 'dx(dy(u))' and 'dy(dx(u))' - atomic_exprs = {SymbolicExpr(a).name : a for a in atomic_exprs} - atomic_exprs = tuple(atomic_exprs.values()) + # Make sure that we only pick one between 'dx1(dx2(u))' and 'dx2(dx1(u))' + new_atoms = {SymbolicExpr(a).name : a for a in new_atoms} + new_atoms = tuple(new_atoms.values()) - for atom in atomic_exprs: + # Create a list of statements for initialization of the point values, + # for each of the atoms in our (possibly new) list. + inits = [] + for atom in new_atoms: + + # Extract field, compute name of coefficient variable, and get base if is_scalar_field(atom): field = atom.atoms(ScalarField).pop() field_name = 'coeff_' + SymbolicExpr(field).name base = field - elif is_vector_field(atom): field = atom.atoms(IndexedVectorField).pop() field_name = 'coeff_' + SymbolicExpr(field).name base = field.base - else: raise TypeError('atom must be either scalar or vector field') - # ... + # Obtain variable for storing point values of test function test_fun = SymbolicExpr(atom.subs(base, test_function)) - # ... - # ... - # TODO: verify if 'get_index_logical_derivatives' should be used instead - orders = [*get_index_derivatives(atom).values()] + # ... + orders = [*get_index(atom).values()] args = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)] inits += [Assign(test_fun, Mul(*args))] # ... @@ -350,7 +371,7 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad, updates += [AugAssign(val,'+',Mul(*args))] # ... - return inits, updates, map_stmts, atomic_exprs + return inits, updates, map_stmts, new_atoms #============================================================================== # TODO: merge into 'compute_atoms_expr_field' From 84f72ade0d7411af45fa66e37b97994c8f9db799 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 29 Aug 2019 19:10:21 +0200 Subject: [PATCH 11/12] [issue #45] Do not use pycode on SymPDE expressions. --- psydac/api/ast/expr.py | 6 ++++-- psydac/api/ast/fem.py | 4 ++-- psydac/api/ast/glt.py | 6 ++++-- psydac/api/expr.py | 18 +++++++--------- psydac/api/printing/pycode.py | 40 ----------------------------------- 5 files changed, 18 insertions(+), 56 deletions(-) diff --git a/psydac/api/ast/expr.py b/psydac/api/ast/expr.py index a8c913a04..a9b9f5d91 100644 --- a/psydac/api/ast/expr.py +++ b/psydac/api/ast/expr.py @@ -340,10 +340,12 @@ def _initialize(self): val = val[indices] if isinstance(expr, (Matrix, ImmutableDenseMatrix)): - body += [Assign(val, expr[i_row,i_col])] + rhs = SymbolicExpr(expr[i_row,i_col]) + body += [Assign(val, rhs)] else: - body += [Assign(val, expr)] + rhs = SymbolicExpr(expr) + body += [Assign(val, rhs)] for i in range(dim-1,-1,-1): x = indices[i] diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index f7443b3e6..56540ac5a 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -678,7 +678,7 @@ def _initialize(self): fields_coeffs = variables(['coeff_{}'.format(f) for f in field_atoms], dtype='real', rank=dim, cls=IndexedVariable) - fields_val = variables(['{}_values'.format(f) for f in fields_str], + fields_val = variables(['{}_values'.format(f) for f in fields_logical_str], dtype='real', rank=dim, cls=IndexedVariable) fields_tmp_coeffs = variables(['tmp_coeff_{}'.format(f) for f in field_atoms], dtype='real', rank=dim, cls=IndexedVariable) @@ -923,7 +923,7 @@ def _initialize(self): for i in range(start, end): if not( i in zero_terms ): - e = Mul(expr[i],wvol) + e = SymbolicExpr(Mul(expr[i],wvol)) body.append(AugAssign(v[i],'+', e)) # ... diff --git a/psydac/api/ast/glt.py b/psydac/api/ast/glt.py index f24251c5c..74da6ec5f 100644 --- a/psydac/api/ast/glt.py +++ b/psydac/api/ast/glt.py @@ -430,10 +430,12 @@ def _initialize(self): symbol = symbol[indices] if isinstance(expr, (Matrix, ImmutableDenseMatrix)): - body += [Assign(symbol, expr[i_row,i_col])] + rhs = SymbolicExpr(expr[i_row,i_col]) + body += [Assign(symbol, rhs)] else: - body += [Assign(symbol, expr)] + rhs = SymbolicExpr(expr) + body += [Assign(symbol, rhs)] for i in range(dim-1,-1,-1): x = indices[i] diff --git a/psydac/api/expr.py b/psydac/api/expr.py index a70192fc0..932e794dd 100644 --- a/psydac/api/expr.py +++ b/psydac/api/expr.py @@ -1,25 +1,23 @@ # coding: utf-8 # TODO for the moment we assume Product of same space +# TODO properly treat expression with mapping -import numpy as np from itertools import product +from sympy import Expr +import numpy as np + +from sympde.expr import TerminalExpr from psydac.api.basic import BasicCodeGen from psydac.api.settings import PSYDAC_BACKEND_PYTHON, PSYDAC_DEFAULT_FOLDER from psydac.api.grid import CollocationBasisValues - from psydac.api.ast.expr import ExprKernel, ExprInterface - from psydac.cad.geometry import Geometry from psydac.mapping.discrete import SplineMapping, NurbsMapping - -from psydac.fem.splines import SplineSpace -from psydac.fem.tensor import TensorFemSpace -from psydac.fem.vector import ProductFemSpace - -from sympde.expr import TerminalExpr -from sympy import Expr +from psydac.fem.splines import SplineSpace +from psydac.fem.tensor import TensorFemSpace +from psydac.fem.vector import ProductFemSpace #============================================================================== class DiscreteExpr(BasicCodeGen): diff --git a/psydac/api/printing/pycode.py b/psydac/api/printing/pycode.py index 685772010..a00314d0f 100644 --- a/psydac/api/printing/pycode.py +++ b/psydac/api/printing/pycode.py @@ -56,45 +56,6 @@ def _print_GltInterface(self, expr): code = '\n'.join(self._print(i) for i in expr.imports) return code +'\n' + self._print(expr.func) - # ......................................................... - - # ......................................................... - # SYMPDE objects - # ......................................................... - def _print_dx(self, expr): - arg = expr.args[0] - if isinstance(arg, _partial_derivatives): - arg = SymbolicExpr(arg).name - else: - arg = self._print(arg) + '_' - return arg + 'x' - - def _print_dy(self, expr): - arg = expr.args[0] - if isinstance(arg, _partial_derivatives): - arg = SymbolicExpr(arg).name - else: - arg = self._print(arg) + '_' - return arg + 'y' - - def _print_dz(self, expr): - arg = expr.args[0] - if isinstance(arg, _partial_derivatives): - arg = SymbolicExpr(arg).name - else: - arg = self._print(arg) + '_' - return arg + 'z' - - def _print_IndexedTestTrial(self, expr): - base = self._print(expr.base) - index = self._print(expr.indices[0]) - return '{base}_{i}'.format(base=base, i=index) - - def _print_IndexedVectorField(self, expr): - base = self._print(expr.base) - index = self._print(expr.indices[0]) - return '{base}_{i}'.format(base=base, i=index) - # ......................................................... # ......................................................... # SYMPY objects @@ -106,7 +67,6 @@ def _print_AppliedUndef(self, expr): args = ','.join(self._print(i) for i in expr.args) fname = self._print(expr.func.__name__) return '{fname}({args})'.format(fname=fname, args=args) - # ......................................................... #============================================================================== def pycode(expr, **settings): From f5be1de9d8bb56fb298f3a495be243cd8bd18b97 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 29 Aug 2019 20:28:30 +0200 Subject: [PATCH 12/12] [issue #45] Add unit test for biharmonic eqn. with mapping. --- .../tests/test_api_1d_compatible_spaces.py | 5 +- .../api/tests/test_api_2d_scalar_mapping.py | 87 +++++++++++++++++++ psydac/api/tests/test_api_expr_2d_scalar.py | 2 +- 3 files changed, 89 insertions(+), 5 deletions(-) diff --git a/psydac/api/tests/test_api_1d_compatible_spaces.py b/psydac/api/tests/test_api_1d_compatible_spaces.py index 14d1f62fe..9c8d83ab2 100644 --- a/psydac/api/tests/test_api_1d_compatible_spaces.py +++ b/psydac/api/tests/test_api_1d_compatible_spaces.py @@ -24,7 +24,6 @@ from psydac.api.discretization import discretize from numpy import linspace, zeros, allclose -import numpy as np from mpi4py import MPI import pytest @@ -97,9 +96,7 @@ def run_system_1_1d_dir(f0, sol, ncells, degree): def test_api_system_1_1d_dir_1(): from sympy.abc import x - + f0 = -(2*pi)**2*sin(2*pi*x) u = sin(2*pi*x) x = run_system_1_1d_dir(f0, u,ncells=[10], degree=[2]) - - diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index 3bd104aaf..6c73888b5 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -264,6 +264,72 @@ def run_laplace_2d_neu(filename, solution, f, comm=None): return l2_error, h1_error +#============================================================================== +def run_biharmonic_2d_dir(filename, solution, f, comm=None): + + # ... abstract model + domain = Domain.from_file(filename) + + V = ScalarFunctionSpace('V', domain) + + F = element_of(V, name='F') + + v = element_of(V, name='v') + u = element_of(V, name='u') + + int_0 = lambda expr: integral(domain , expr) + + expr = laplace(v) * laplace(u) + a = BilinearForm((v,u),int_0(expr)) + + expr = f*v + l = LinearForm(v, int_0(expr)) + + error = F - solution + l2norm = Norm(error, domain, kind='l2') + h1norm = Norm(error, domain, kind='h1') + h2norm = Norm(error, domain, kind='h2') + + nn = NormalVector('nn') + bc = [EssentialBC(u, 0, domain.boundary)] + bc += [EssentialBC(dot(grad(u), nn), 0, domain.boundary)] + equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) + # ... + + # ... create the computational domain from a topological domain + domain_h = discretize(domain, filename=filename, comm=comm) + # ... + + # ... discrete spaces + Vh = discretize(V, domain_h) + # ... + + # ... dsicretize the equation using Dirichlet bc + equation_h = discretize(equation, domain_h, [Vh, Vh]) + # ... + + # ... discretize norms + l2norm_h = discretize(l2norm, domain_h, Vh) + h1norm_h = discretize(h1norm, domain_h, Vh) + h2norm_h = discretize(h2norm, domain_h, Vh) + # ... + + # ... solve the discrete equation + x = equation_h.solve() + # ... + + # ... + phi = FemField( Vh, x ) + # ... + + # ... compute norms + l2_error = l2norm_h.assemble(F=phi) + h1_error = h1norm_h.assemble(F=phi) + h2_error = h2norm_h.assemble(F=phi) + # ... + + return l2_error, h1_error, h2_error + ############################################################################### # SERIAL TESTS ############################################################################### @@ -521,6 +587,27 @@ def test_api_laplace_2d_neu_identity(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) +#============================================================================== +def test_api_biharmonic_2d_dir_identity(): + filename = os.path.join(mesh_dir, 'quart_circle.h5') + + from sympy.abc import x,y + + solution = (sin(pi*x)*sin(pi*y))**2 + f = laplace(laplace(solution)) + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f) + + expected_l2_error = 0.3092792236008727 + expected_h1_error = 1.3320441589030887 + expected_h2_error = 6.826223014197719 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + assert( abs(h2_error - expected_h2_error) < 1.e-7) + +# TODO: add biharmonic on Collela and quart_circle mappings + ##============================================================================== ## TODO DEBUG, not working since merge with devel #def test_api_laplace_2d_neu_collela(): diff --git a/psydac/api/tests/test_api_expr_2d_scalar.py b/psydac/api/tests/test_api_expr_2d_scalar.py index d6aef4024..b88da5720 100644 --- a/psydac/api/tests/test_api_expr_2d_scalar.py +++ b/psydac/api/tests/test_api_expr_2d_scalar.py @@ -69,6 +69,6 @@ def run_poisson_2d_dir(ncells, degree, comm=None): ############################################################################### #============================================================================== -def test_api_glt_poisson_2d_dir_1(): +def test_api_expr_2d_1(): run_poisson_2d_dir(ncells=[2**3,2**3], degree=[2,2])