diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py index f7855277d79e9ff5b94d2effde9a98a7df597776..4ce2f22a8c640aca9f5a6b60acdd497182493bd2 100644 --- a/python/dune/perftool/compile.py +++ b/python/dune/perftool/compile.py @@ -9,6 +9,7 @@ def read_ufl(uflfile): from ufl.algorithms.formfiles import read_ufl_file, interpret_ufl_namespace from ufl.algorithms import compute_form_data + # Read the given ufl file and execute it uflcode = read_ufl_file(uflfile) namespace = {} uflcode = "from dune.perftool.ufl.execution import *\n" + uflcode @@ -23,23 +24,33 @@ def read_ufl(uflfile): with file(pyname, "w") as f: f.write(pycode) raise SyntaxError("Not a valid ufl file, dumped a debug script: {}".format(pyname)) + + # Extract the form from the given data data = interpret_ufl_namespace(namespace) - formdata = compute_form_data(data.forms[0]) + form = data.forms[0] + form = compute_form_data(form).preprocessed_form # We do not expect more than one form assert len(data.forms) == 1 - # We make some assumptions on the UFL expression! + # Before we go any deeper, try check validity (although it does basically the same things we do later) + # TODO: do not make this an assertion of course! from dune.perftool.ufl.validity import check_validity - check_validity(data.forms[0]) + assert check_validity(form) + + # apply some transformations unconditionally! + from dune.perftool.ufl.transformations import transform_form + from dune.perftool.ufl.transformations.splitarguments import split_arguments + + form = transform_form(form, split_arguments) - return formdata + return form -def generate_driver(formdata, filename): +def generate_driver(form, filename): # The driver module uses a global variable for ease of use - from dune.perftool.pdelab.driver import set_formdata - set_formdata(formdata) + from dune.perftool.pdelab.driver import set_form + set_form(form) # The vtkoutput is the generating method that triggers all others. # Alternatively, one could use the `solve` method. @@ -73,11 +84,11 @@ def generate_driver(formdata, filename): def compile_form(): from dune.perftool.options import get_option - ufldata = read_ufl(get_option("uflfile")) + form = read_ufl(get_option("uflfile")) if get_option("driver_file"): - generate_driver(ufldata, get_option("driver_file")) + generate_driver(form, get_option("driver_file")) if get_option("operator_file"): from dune.perftool.pdelab.localoperator import generate_localoperator - generate_localoperator(ufldata, get_option("operator_file")) + generate_localoperator(form, get_option("operator_file")) diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py new file mode 100644 index 0000000000000000000000000000000000000000..d7c5ed53ebdd4fc0a528cbf16241151e7f197d0a --- /dev/null +++ b/python/dune/perftool/loopy/transformer.py @@ -0,0 +1,160 @@ +""" +This is the module that contains the main transformation from ufl to loopy +(with pdelab as the hardcoded generation target) +""" + +from __future__ import absolute_import +from ufl.algorithms import MultiFunction +# Spread the pymbolic import statements to where they are used. +from pymbolic.primitives import Variable, Subscript, Sum, Product +import loopy +import numpy +import ufl + +from dune.perftool import Restriction + +# Define the generators that are used here +from dune.perftool.generation import generator_factory +loopy_iname = generator_factory(item_tags=("loopy", "kernel", "iname")) +loopy_expr_instruction = generator_factory(item_tags=("loopy", "kernel", "instruction", "exprinstruction"), no_deco=True) +loopy_temporary_variable = generator_factory(item_tags=("loopy", "kernel", "temporary"), on_store=lambda n: loopy.TemporaryVariable(n, dtype=numpy.float64), no_deco=True) +loopy_c_instruction = generator_factory(item_tags=("loopy", "kernel", "instruction", "cinstruction"), no_deco=True) +loopy_valuearg = generator_factory(item_tags=("loopy", "kernel", "argument", "valuearg"), on_store=lambda n: loopy.ValueArg(n), no_deco=True) + +@generator_factory(item_tags=("loopy", "kernel", "argument", "globalarg")) +def loopy_globalarg(name, shape=loopy.auto): + if isinstance(shape, str): + shape = (shape,) + return loopy.GlobalArg(name, numpy.float64, shape) + +@generator_factory(item_tags=("loopy", "kernel", "domain")) +def loopy_domain(iname, shape): + loopy_valuearg(shape) + return "{{ [{0}] : 0<={0}<{1} }}".format(iname, shape) + +@loopy_iname +def dimension_iname(index): + from dune.perftool.pdelab import name_index + from dune.perftool.pdelab.geometry import name_dimension + iname = name_index(index) + dimname = name_dimension() + loopy_domain(iname, dimname) + return iname + +@loopy_iname +def argument_iname(arg): + # TODO extract the {iname}_n thing by a preamble + from dune.perftool.ufl.modified_terminals import ModifiedArgumentNumber + iname = "arg{}".format(chr(ord("i") + ModifiedArgumentNumber()(arg))) + loopy_domain(iname, iname + "_n") + return iname + +@loopy_iname +def quadrature_iname(): + loopy_domain("q", "q_n") + return "q" + +from dune.perftool.ufl.modified_terminals import ModifiedTerminalTracker +from dune.perftool.pymbolic.uflmapper import UFL2PymbolicMapper + +class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper): + def __init__(self, embracing_mf): + super(UFL2LoopyVisitor, self).__init__() + self.call = embracing_mf.call + + def coefficient(self, o): + # All trial functions should already be handled by the TrialFunctionExtractor class! + assert o.count() != 0 + # TODO implement non-trialfunction coefficients + raise NotImplementedError + + # TODO use multiple inheritance and have a geometry transformer in the pdelab subpackage + def facet_area(self, o): + from pymbolic.primitives import Variable + from dune.perftool.pdelab.geometry import name_facetarea + return Variable(name_facetarea()) + + +# TODO This pattern here could be abstracted as some sort of "SupersedingMultiFunction": +# A superseding multifunction is always asked to handle the given expression first. Only +# if it explicitly delegates, the underlying mf is called. Whenever, the underlying mf +# goes into recursion, it asks back to the superseding function. Martin would pity me for +# the performance overhead, but it makes life easier. +class TrialFunctionExtractor(MultiFunction): + call = MultiFunction.__call__ + + def __call__(self, o): + from dune.perftool.ufl.modified_terminals import ModifiedArgumentExtractor + self.tf = ModifiedArgumentExtractor()(o, trialfunction=True) + self.u2l = UFL2LoopyVisitor(self) + + return self.call(o) + + def expr(self, o): + if o in self.tf: + # This is a modified trial function! + from dune.perftool.pdelab.argument import name_trialfunction + name = name_trialfunction(o) + loopy_globalarg(name) + return Variable(name) + else: + return self.u2l(o) + +class _Counter: + counter = 0 + +def get_count(): + c = _Counter.counter + _Counter.counter = c + 1 + return c + +def transform_accumulation_term(term): + # Get the accumulation expression and the modified arguments + expr, args = term + + # We always have a quadrature loop + quadrature_iname() + + # Get the pymbolic expression needed for this accumulation term. + # This includes filling the cache with all sorts of necessary preambles! + pymbolic_expr = TrialFunctionExtractor()(expr) + + # Now simplify the expression + # TODO: Add a switch to disable/configure this. +# from dune.perftool.pymbolic.simplify import simplify_pymbolic_expression +# pymbolic_expr = simplify_pymbolic_expression(pymbolic_expr) + + # Define a temporary variable for this expression + expr_tv_name = "expr_"+str(get_count()).zfill(4) + expr_tv = loopy_temporary_variable(expr_tv_name) + loopy_expr_instruction(loopy.ExpressionInstruction(assignee=Variable(expr_tv_name), expression=pymbolic_expr)) + + # The data that is used to collect the arguments for the accumulate function + accumargs = [] + argument_code = [] + + # Generate the code for the modified arguments: + for arg in args: + from dune.perftool.pdelab.argument import name_argumentspace, name_argument + accumargs.append(name_argumentspace(arg)) + accumargs.append(argument_iname(arg)) + name = name_argument(arg) + argument_code.append(name) + loopy_globalarg(name) + + from dune.perftool.pdelab.argument import name_residual + residual = name_residual() + + from dune.perftool.generation import retrieve_cache_items + inames = retrieve_cache_items("iname") + + from dune.perftool.pdelab.quadrature import name_factor + loopy_c_instruction(loopy.CInstruction(inames, + "{}.accumulate({}, {}*{}*{})".format(residual, + ", ".join(accumargs), + expr_tv_name, + "*".join(argument_code), + name_factor() + ) + ) + ) diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py index 1e0d1e89e6538e59e801455357fde30041d05a0c..8809c1a5fa1d3f2c1808991d1c0e8c457250c930 100644 --- a/python/dune/perftool/options.py +++ b/python/dune/perftool/options.py @@ -27,4 +27,8 @@ def get_form_compiler_arguments(): return args def get_option(key, default=None): - return get_form_compiler_arguments().get(key, default) \ No newline at end of file + try: + __IPYTHON__ + return default + except NameError: + return get_form_compiler_arguments().get(key, default) \ No newline at end of file diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py index 23c40d5a4c28155c32d02d8ba0a6d74c822f9efb..2d98c18e408bf09684ed404bb1a5b7d9f81bf604 100644 --- a/python/dune/perftool/pdelab/__init__.py +++ b/python/dune/perftool/pdelab/__init__.py @@ -6,7 +6,7 @@ dune_symbol = generator_factory(item_tags=("pdelab", "kernel", "symbol")) dune_preamble = generator_factory(item_tags=("pdelab", "kernel", "preamble"), counted=True) dune_include = generator_factory(on_store=lambda i: "#include<{}>".format(i), item_tags=("pdelab", "include"), no_deco=True) -from dune.perftool.transformer import quadrature_iname +from dune.perftool.loopy.transformer import quadrature_iname from loopy import CInstruction def quadrature_preamble(assignees=[]): diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py index 0e951fed920aa61ea5e311ebbf5c6efd32f0eaaa..71d989cb85c18e9d870da1a25a8393c0d59f37b1 100644 --- a/python/dune/perftool/pdelab/argument.py +++ b/python/dune/perftool/pdelab/argument.py @@ -1,17 +1,19 @@ """ Generator functions related to trial and test functions and the accumulation loop""" from dune.perftool.pdelab import dune_symbol, dune_preamble - +from dune.perftool.ufl.modified_terminals import ModifiedArgumentDescriptor @dune_symbol -def name_testfunction(arg, index, grad, restriction): - if len(arg.ufl_element().sub_elements()) > 0: +def name_testfunction(modarg): + ma = ModifiedArgumentDescriptor(modarg) + if len(ma.expr.ufl_element().sub_elements()) > 0: pass - return "{}a{}".format("grad_" if grad else "", arg.number()) + return "{}a{}".format("grad_" if ma.grad else "", ma.expr.number()) @dune_symbol -def name_trialfunction(arg, index, grad, restriction): - return "{}c{}".format("grad_" if grad else "", arg.count()) +def name_trialfunction(modarg): + ma = ModifiedArgumentDescriptor(modarg) + return "{}c{}".format("grad_" if ma.grad else "", ma.expr.count()) @dune_symbol def name_testfunctionspace(*a): @@ -23,11 +25,22 @@ def name_trialfunctionspace(*a): #TODO return "lfsu" -def name_argumentspace(arg, *others): - if arg.number() == 0: - return name_testfunctionspace(arg, *others) - if arg.number() == 1: - return name_trialfunctionspace(arg, *others) +def name_argumentspace(modarg): + ma = ModifiedArgumentDescriptor(modarg) + if ma.expr.number() == 0: + return name_testfunctionspace(modarg) + if ma.expr.number() == 1: + return name_trialfunctionspace(modarg) + # We should never encounter an argument other than 0 or 1 + assert False + + +def name_argument(modarg): + ma = ModifiedArgumentDescriptor(modarg) + if ma.expr.number() == 0: + return name_testfunction(modarg) + if ma.expr.number() == 1: + return name_trialfunction(modarg) # We should never encounter an argument other than 0 or 1 assert False diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py index 92f3549454f1dede277d7cc81cb28f954551df5b..4dd55ca8d95e68f4aef951ed7e52a164f1d07a25 100644 --- a/python/dune/perftool/pdelab/driver.py +++ b/python/dune/perftool/pdelab/driver.py @@ -12,12 +12,12 @@ from dune.perftool.pdelab import dune_symbol, dune_preamble, dune_include # Have a global variable with the entire form data. This allows functions that depend # deterministically on the entire data set to directly access it instead of passing it # through the entire generator chain. -_formdata = None +_form = None # Have a function access this global data structure -def set_formdata(fd): - global _formdata - _formdata = fd +def set_form(f): + global _form + _form = f def isLagrange(fem): @@ -73,7 +73,7 @@ def name_initree(): @dune_preamble def define_dimension(name): - return "static const int {} = {};".format(name, _formdata.geometric_dimension) + return "static const int {} = {};".format(name, __form.geometric_dimension) @dune_symbol def name_dimension(): @@ -83,11 +83,11 @@ def name_dimension(): @dune_preamble def typedef_grid(name): dim = name_dimension() - if any(_formdata.unique_elements[0].cell().cellname() in x for x in ["vertex", "interval", "quadrilateral", "hexalateral"]): + if any(_form.unique_elements[0].cell().cellname() in x for x in ["vertex", "interval", "quadrilateral", "hexalateral"]): gridt = "Dune::YaspGrid<{}>".format(dim) dune_include("dune/grid/yaspgrid.hh") else: - if any(_formdata.unique_elements[0].cell().cellname() in x for x in ["triangle", "tetrahedron"]): + if any(_form.unique_elements[0].cell().cellname() in x for x in ["triangle", "tetrahedron"]): gridt = "Dune::UGGrid<{}>".format(dim) dune_include("dune/grid/uggrid.hh") else: @@ -323,11 +323,11 @@ def name_gfs(expr): @dune_preamble def define_dofestimate(name): # Provide a worstcase estimate for the number of entries per row based on the given gridfunction space and cell geometry - if isQuadrilateral(_formdata.coefficient_elements[0]): + if isQuadrilateral(_form.coefficient_elements[0]): geo_factor = "4" else: geo_factor = "6" - gfs = name_gfs(_formdata.coefficient_elements[0]) + gfs = name_gfs(_form.coefficient_elements[0]) ini = name_initree() return ["int generic_dof_estimate = {} * {}.maxLocalSize();".format(geo_factor, gfs), "int dof_estimate = {}.get<int>(\"istl.number_of_nnz\", generic_dof_estimate);".format(ini)] @@ -401,11 +401,11 @@ def name_localoperator(): @dune_preamble def typedef_gridoperator(name): - ugfs = type_gfs(_formdata.coefficient_elements[0]) - vgfs = type_gfs(_formdata.argument_elements[0]) + ugfs = type_gfs(_form.coefficient_elements[0]) + vgfs = type_gfs(_form.argument_elements[0]) lop = type_localoperator() - ucc = type_constraintscontainer(_formdata.coefficient_elements[0]) - vcc = type_constraintscontainer(_formdata.argument_elements[0]) + ucc = type_constraintscontainer(_form.coefficient_elements[0]) + vcc = type_constraintscontainer(_form.argument_elements[0]) mb = type_matrixbackend() df = type_domainfield() r = type_range() @@ -420,10 +420,10 @@ def type_gridoperator(): @dune_preamble def define_gridoperator(name): gotype = type_gridoperator() - ugfs = name_gfs(_formdata.coefficient_elements[0]) - ucc = name_constraintscontainer(_formdata.coefficient_elements[0]) - vgfs = name_gfs(_formdata.argument_elements[0]) - vcc = name_constraintscontainer(_formdata.argument_elements[0]) + ugfs = name_gfs(_form.coefficient_elements[0]) + ucc = name_constraintscontainer(_form.coefficient_elements[0]) + vgfs = name_gfs(_form.argument_elements[0]) + vcc = name_constraintscontainer(_form.argument_elements[0]) lop = name_localoperator() mb = name_matrixbackend() return "{} {}({}, {}, {}, {}, {}, {});".format(gotype, name, ugfs, ucc, vgfs, vcc, lop, mb) @@ -446,7 +446,7 @@ def type_vector(): @dune_preamble def define_vector(name): vtype = type_vector() - gfs = name_gfs(_formdata.coefficient_elements[0]) + gfs = name_gfs(_form.coefficient_elements[0]) return ["{} {}({});".format(vtype, name, gfs), "{} = 0.0;".format(name)] @dune_symbol @@ -514,7 +514,7 @@ def name_stationarylinearproblemsolver(): @dune_preamble def dune_solve(): from ufl.algorithms.predicates import is_multilinear - if is_multilinear(_formdata.preprocessed_form): + if is_multilinear(_form.preprocessed_form): slp = name_stationarylinearproblemsolver() return "{}.apply();".format(slp) else: @@ -535,7 +535,7 @@ def name_vtkfile(): def vtkoutput(): dune_include("dune/pdelab/gridfunctionspace/vtk.hh") vtkwriter = name_vtkwriter() - gfs = name_gfs(_formdata.coefficient_elements[0]) + gfs = name_gfs(_form.coefficient_elements[0]) vec = name_vector() vtkfile = name_vtkfile() dune_solve() diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py index bc2bc74c37d6dec90897efba3d8e42e63ea93ce1..6d0db436c7cae4aed6d9d021f5097c8ee8f3aa36 100644 --- a/python/dune/perftool/pdelab/geometry.py +++ b/python/dune/perftool/pdelab/geometry.py @@ -4,3 +4,8 @@ from dune.perftool.pdelab import dune_symbol def name_dimension(): #TODO preamble define_dimension return "dim" + +@dune_symbol +def name_facetarea(): + #TODO preambles + return "farea" \ No newline at end of file diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py index 2cada9447ae045d5db7deeec1239363b9124e94b..c7885b87fd0b3d9bb141a6664d0f0d57f3840666 100644 --- a/python/dune/perftool/pdelab/localoperator.py +++ b/python/dune/perftool/pdelab/localoperator.py @@ -1,3 +1,4 @@ +from __future__ import absolute_import from pytools import memoize from dune.perftool.options import get_option from dune.perftool.generation import generator_factory @@ -78,9 +79,23 @@ def generate_term(integrand=None, measure=None): # Get the measure specifics specifics = measure_specific_details(measure) - # Transform the expression. All relevant results are then in the cache - from dune.perftool.transformer import transform_expression - transform_expression(integrand) + from dune.perftool.ufl.rank import ufl_rank + if ufl_rank(integrand) == 2: + from IPython import embed; embed() + + # Now split the given integrand into accumulation expressions + from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms + accterms = split_into_accumulation_terms(integrand) + + # Iterate over the terms and generate a kernel + for term in accterms: + from dune.perftool.loopy.transformer import transform_accumulation_term + transform_accumulation_term(term) + +# +# # Transform the expression. All relevant results are then in the cache +# from dune.perftool.loopy.transformer import transform_expression +# transform_expression(integrand) # Extract the information, which is needed to create a loopy kernel. # First extracting it, might be useful to alter it before kernel generation. @@ -102,11 +117,7 @@ def generate_term(integrand=None, measure=None): return str(generate_code(kernel)[0]) -def generate_localoperator(ufldata, operatorfile): - form = ufldata.preprocessed_form - - operator_methods = [] - +def generate_localoperator(form, operatorfile): # For the moment, I do assume that there is but one integral of each type. This might differ # if you use different quadrature orders for different terms. assert len(form.integrals()) == len(set(i.integral_type() for i in form.integrals())) @@ -115,6 +126,9 @@ def generate_localoperator(ufldata, operatorfile): from dune.perftool.generation import delete_cache delete_cache() + # Have a data structure collect the generated kernels + operator_methods = [] + # Generate the necessary residual methods for integral in form.integrals(): body = generate_term(integrand=integral.integrand(), measure=integral.integral_type()) diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py index 06b167541b98b9e9baa089c2fc133febcbed0fae..e306e3943e6581d7cfd87c0a5110657a10c7a2a7 100644 --- a/python/dune/perftool/pdelab/quadrature.py +++ b/python/dune/perftool/pdelab/quadrature.py @@ -1,5 +1,5 @@ from dune.perftool.generation import generator_factory -from dune.perftool.transformer import quadrature_iname, loopy_temporary_variable +from dune.perftool.loopy.transformer import quadrature_iname, loopy_temporary_variable from dune.perftool.pdelab import dune_symbol, quadrature_preamble, dune_preamble @dune_symbol diff --git a/python/dune/perftool/pymbolic/__init__.py b/python/dune/perftool/pymbolic/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/python/dune/perftool/pymbolic/simplify.py b/python/dune/perftool/pymbolic/simplify.py new file mode 100644 index 0000000000000000000000000000000000000000..3382cc245540971ceb00315244e2bf5da90edcf3 --- /dev/null +++ b/python/dune/perftool/pymbolic/simplify.py @@ -0,0 +1,9 @@ +from __future__ import absolute_import + +from pymbolic.sympy_interface import SympyToPymbolicMapper, PymbolicToSympyMapper +from sympy import simplify + +def simplify_pymbolic_expression(e): + sympyexpr = PymbolicToSympyMapper()(e) + simplified = simplify(sympyexpr) + return SympyToPymbolicMapper()(simplified) \ No newline at end of file diff --git a/python/dune/perftool/pymbolic/uflmapper.py b/python/dune/perftool/pymbolic/uflmapper.py new file mode 100644 index 0000000000000000000000000000000000000000..535c2af747574ecb675a4b665e193c85ba84bde8 --- /dev/null +++ b/python/dune/perftool/pymbolic/uflmapper.py @@ -0,0 +1,59 @@ +""" A multi function mapping UFL expressions to pymbolic expressions + +This mapper only defines handler for those UFL expression types that have an +equivalent in pymbolic (mainly algebraic operations). + +This is mainly intended as a base class for anything mapping ufl to pymbolic. +""" + +from __future__ import absolute_import +from ufl.algorithms import MultiFunction +from pymbolic.primitives import Product, Quotient, Subscript, Sum, Variable + +# The constructed pymbolic expressions use n-ary operators instead of ufls binary operators +from dune.perftool.ufl.flatoperators import get_operands + +class UFL2PymbolicMapper(MultiFunction): + def __init__(self): + super(UFL2PymbolicMapper, self).__init__() + + call = MultiFunction.__call__ + + def product(self, o): + return Product(tuple(self.call(op) for op in get_operands(o))) + + def multi_index(self, o): + return tuple(self.call(op) for op in o.ufl_operands) + + def index(self, o): + # One might as well take the uflname as string here, but I apply this function + from dune.perftool.pdelab import name_index + + return Variable(name_index(o)) + + def fixed_index(self, o): + return o._value + + def indexed(self, o): + aggr = self.call(o.ufl_operands[0]) + ind = self.call(o.ufl_operands[1]) + + # If the left operand is already a subscripted pymbolic expression, we need to do tuple addition on the indices + if isinstance(aggr, Subscript): + return Subscript(aggr.aggregate, aggr.index + ind) + else: + return Subscript(aggr, ind) + + def float_value(self, o): + return o.value() + + def int_value(self, o): + return o.value() + + def division(self, o): + assert len(o.ufl_operands) == 2 + + return Quotient(self.call(o.ufl_operands[0]), self.call(o.ufl_operands[1])) + + def sum(self, o): + return Sum(tuple(self.call(op) for op in get_operands(o))) diff --git a/python/dune/perftool/transformer.py b/python/dune/perftool/transformer.py deleted file mode 100644 index 55b042d70cad63c2b2bb19282b0916dac8de87b3..0000000000000000000000000000000000000000 --- a/python/dune/perftool/transformer.py +++ /dev/null @@ -1,194 +0,0 @@ -""" -This is the module that contains the main transformation from ufl to loopy -(with pdelab as the hardcoded generation target) -""" - -from __future__ import absolute_import -from ufl.algorithms import MultiFunction -# Spread the pymbolic import statements to where they are used. -from pymbolic.primitives import Variable, Subscript, Sum, Product -import loopy -import numpy -import ufl - -from dune.perftool import Restriction - -# Define the generators that are used here -from dune.perftool.generation import generator_factory -loopy_iname = generator_factory(item_tags=("loopy", "kernel", "iname")) -loopy_expr_instruction = generator_factory(item_tags=("loopy", "kernel", "instruction", "exprinstruction"), no_deco=True) -loopy_temporary_variable = generator_factory(item_tags=("loopy", "kernel", "temporary"), on_store=lambda n: loopy.TemporaryVariable(n, dtype=numpy.float64), no_deco=True) -loopy_c_instruction = generator_factory(item_tags=("loopy", "kernel", "instruction", "cinstruction"), no_deco=True) -loopy_valuearg = generator_factory(item_tags=("loopy", "kernel", "argument", "valuearg"), on_store=lambda n: loopy.ValueArg(n), no_deco=True) - -@generator_factory(item_tags=("loopy", "kernel", "argument", "globalarg")) -def loopy_globalarg(name, shape=loopy.auto): - if isinstance(shape, str): - shape = (shape,) - return loopy.GlobalArg(name, numpy.float64, shape) - -@generator_factory(item_tags=("loopy", "kernel", "domain")) -def loopy_domain(iname, shape): - loopy_valuearg(shape) - return "{{ [{0}] : 0<={0}<{1} }}".format(iname, shape) - -@loopy_iname -def dimension_iname(index): - from dune.perftool.pdelab import name_index - from dune.perftool.pdelab.geometry import name_dimension - iname = name_index(index) - dimname = name_dimension() - loopy_domain(iname, dimname) - return iname - -@loopy_iname -def argument_iname(arg): - # TODO extract the {iname}_n thing by a preamble - iname = "arg{}".format(chr(ord("i") + arg.number())) - loopy_domain(iname, iname + "_n") - return iname - -@loopy_iname -def quadrature_iname(): - loopy_domain("q", "q_n") - return "q" - - -class UFLVisitor(MultiFunction): - def __call__(self, o, init_inames): - # Make this multifunction stateful: Store information similar to uflacs' modified terminal - self.grad = False - self.reference_grad = False - self.index = None - self.restriction = Restriction.NONE - - # Have a short cut for the bases class call operator - self.call = lambda *a : MultiFunction.__call__(self, *a) - - # Have a list of arguments that this term depends on. - self.arguments = [] - self.inames = init_inames - - # Visit the expression. The return value will be a pymbolic expression. - loopyexpr = self.call(o) - - # The expression is completely processes, an accumulation instruction can be issued! - # Have some data structures for the analysis of the given arguments. - accumargs = [] - - # visit the given arguments and extract necessary information - for arg, index, grad, restriction in sorted(self.arguments, key=lambda x: x[0].number()): - from dune.perftool.pdelab.argument import name_argumentspace - accumargs.append(name_argumentspace(arg, index, restriction)) - accumargs.append(argument_iname(arg)) - - # Explicitly state the dependency on the quadrature iname - self.inames.append(quadrature_iname()) - - # Evaluate the update factor - from dune.perftool.pdelab.quadrature import name_factor - loopyexpr = Product((loopyexpr, Variable(name_factor()))) - # TODO: Some unique name mangling on this temporary variable is necessary! - expr_tv = loopy_temporary_variable("expr") - loopy_expr_instruction(loopy.ExpressionInstruction(assignee=Variable("expr"), expression=loopyexpr)) - - from dune.perftool.pdelab.argument import name_residual - residual = name_residual() - - loopy_c_instruction(loopy.CInstruction(self.inames, "{}.accumulate({}, {})".format(residual, ", ".join(accumargs), Variable("xyz")))) - - - def grad(self, o): - assert len(o.operands()) == 1 - - self.grad = True - ret = self.call(o.operands()[0]) - self.grad = False - return ret - - def argument(self, o): - assert len(o.operands()) == 0 - - # Construct a tuple that represents this argument and all its modifications - self.arguments.append((o, self.index, self.grad, self.restriction)) - - # Generate the variable name for the test function - from dune.perftool.pdelab.argument import name_testfunction - name = name_testfunction(o, self.index, self.grad, self.restriction) - index = argument_iname(o) - self.inames.append(index) - loopy_globalarg(name) - - return Subscript(Variable(name), (Variable(index),)) - - def coefficient(self, o): - assert len(o.operands()) == 0 - # TODO: Implement coefficient functions - assert o.count() == 0 - - # TODO implement non-trialfunction coefficients! - from dune.perftool.pdelab.argument import name_trialfunction - from dune.perftool.pdelab import name_index - name = name_trialfunction(o, self.index, self.grad, self.restriction) - loopy_globalarg(name) - - return Variable(name) - - def product(self, o): - return Product(tuple(self.call(op) for op in o.operands())) - - def multi_index(self, o): - # I don't think we should ever find a multi index, because its father should take care of it. - raise ValueError - - def indexed(self, o): - # TODO currently, we assume single indices - assert len(o.operands()[1]) == 1 - - # We have two types of indices in our expressions: - # * those that should result in variable subscripts - # * those that should result in function space child selection - # I think it is correct to assume that Indexed expressions - # are of the latter type if and only if their first operand is - # an argument or coefficient. - if isinstance(o.operands()[0], (ufl.Argument, ufl.Coefficient)): - assert not self.index - self.index = o.operands()[1] - ret = self.call(o.operands()[0]) - self.index = None - return ret - else: - from dune.perftool.pdelab import name_index - ret = self.call(o.operands()[0]) - index = Variable(name_index(o.operands()[1])) - if isinstance(ret, Subscript): - return Subscript(ret.aggregate, ret.index + (index,)) - else: - return Subscript(ret, (index,)) - - -class TopSumSeparation(MultiFunction): - """ A multifunction that separates the toplevel sum """ - def __init__(self, visitor=UFLVisitor()): - MultiFunction.__init__(self) - self.visitor = visitor - self.inames = [] - - def expr(self, o): - self.visitor(o, self.inames) - self.inames = [] - - def sum(self, o): - for op in o.operands(): - self(op) - - def index_sum(self, o): - # Generate an iname: We do this here and restrict ourselves to shape dim. - # TODO revisit when implementing general index sums - if isinstance(self.visitor, UFLVisitor): - self.inames.append(dimension_iname(o.operands()[1])) - self(o.operands()[0]) - - -def transform_expression(expr): - TopSumSeparation()(expr) diff --git a/python/dune/perftool/ufl/flatoperators.py b/python/dune/perftool/ufl/flatoperators.py new file mode 100644 index 0000000000000000000000000000000000000000..911d9bde2908c1551c150bfbbb3a8fb02a900e70 --- /dev/null +++ b/python/dune/perftool/ufl/flatoperators.py @@ -0,0 +1,30 @@ +""" A module that treats the problem of binary ufl operators """ + + +def get_operands(expr, operator=None): + """ Get the operators of the flattened operation in the expression """ + # If no operator is specified, it is the root of the expression + if not operator: + operator = type(expr) + + # We assume all operators to be binary + assert len(expr.ufl_operands) == 2 + + result = [] + for op in expr.ufl_operands: + if isinstance(op, operator): + result.extend(get_operands(op, operator)) + else: + result.append(op) + + return result + + +def construct_binary_operator(operands, operator): + if len(operands) == 1: + return operands[0] + if len(operands) == 2: + return operator(operands[0], operands[1]) + else: + mid = len(operands)//2 + 1 + return operator(construct_binary_operator(operands[:mid], operator), construct_binary_operator(operands[mid:len(operands)], operator)) diff --git a/python/dune/perftool/ufl/modified_terminals.py b/python/dune/perftool/ufl/modified_terminals.py index 00eb30da83e555782faebf2257fb895b8a1dabe5..9c7578231294cf5d1bb5ab49e395b20120e02645 100644 --- a/python/dune/perftool/ufl/modified_terminals.py +++ b/python/dune/perftool/ufl/modified_terminals.py @@ -5,6 +5,10 @@ from dune.perftool import Restriction class ModifiedTerminalTracker(MultiFunction): + """ A multifunction base that defines handler for + grad, reference_grad, positive_restricted and negative_restricted. + The appearance of those classes changes the internal state of the MF. + """ def __init__(self): MultiFunction.__init__(self) self.grad = False @@ -12,25 +16,129 @@ class ModifiedTerminalTracker(MultiFunction): self.restriction = Restriction.NONE def positive_restricted(self, o): + assert self.restriction == Restriction.NONE self.restriction = Restriction.POSITIVE - ret = self.expr(o) + ret = self.call(o.ufl_operands[0]) self.restriction = Restriction.NONE return ret def negative_restricted(self, o): + assert self.restriction == Restriction.NONE self.restriction = Restriction.NEGATIVE - ret = self.expr(o) + ret = self.call(o.ufl_operands[0]) self.restriction = Restriction.NONE return ret def grad(self, o): + assert not self.grad and not self.reference_grad self.grad = True - ret = self.expr(o) + ret = self.call(o.ufl_operands[0]) self.grad = False return ret def reference_grad(self, o): + assert not self.reference_grad and not self.grad self.reference_grad = True - ret = self.expr(o) + ret = self.call(o.ufl_operands[0]) self.reference_grad = False return ret + + +class ModifiedArgumentExtractor(MultiFunction): + """ A multifunction that extracts and returns the set of modified arguments """ + + def __call__(self, o, argnumber=None, trialfunction=False): + self.argnumber = argnumber + self.trialfunction = trialfunction + self.modified_arguments = set() + ret = self.call(o) + if ret: + # This indicates that this entire expression was a modified thing... + self.modified_arguments.add(ret) + return tuple(self.modified_arguments) + + def expr(self, o): + for op in o.ufl_operands: + ret = self.call(op) + if ret: + self.modified_arguments.add(ret) + + def indexed(self, o): + if self.call(o.ufl_operands[0]): + return o + + def positive_restricted(self, o): + if self.call(o.ufl_operands[0]): + return o + + def negative_restricted(self, o): + if self.call(o.ufl_operands[0]): + return o + + def grad(self, o): + if self.call(o.ufl_operands[0]): + return o + + def reference_grad(self, o): + if self.call(o.ufl_operands[0]): + return o + + def argument(self, o): + if not self.trialfunction: + if self.argnumber is None or o.number() == self.argnumber: + return o + + def coefficient(self, o): + if self.trialfunction: + if o.count() == 0: + return o + + call = MultiFunction.__call__ + + +class ModifiedArgumentNumber(MultiFunction): + """ return the number() of a modified argument """ + def expr(self, o): + return self(o.ufl_operands[0]) + + def argument(self, o): + return o.number() + + +class ModifiedArgumentDescriptor(MultiFunction): + def __init__(self, e): + MultiFunction.__init__(self) + + self.grad = False + self.reference_grad = False + self.index = None + self.restriction = Restriction.NONE + + self.__call__(e) + self.__call__ = None + + def grad(self, o): + self.grad = True + self(o.ufl_operands[0]) + + def reference_grad(self, o): + self.reference_grad = True + self(o.ufl_operands[0]) + + def positive_restricted(self, o): + self.restriction = Restriction.POSITIVE + self(o.ufl_operands[0]) + + def negative_restricted(self, o): + self.restriction = Restriction.NEGATIVE + self(o.ufl_operands[0]) + + def indexed(self, o): + indexed = o.ufl_operands[1] + self(o.ufl_operands[0]) + + def argument(self, o): + self.expr = o + + def coefficient(self, o): + self.expr = o diff --git a/python/dune/perftool/ufl/rank.py b/python/dune/perftool/ufl/rank.py new file mode 100644 index 0000000000000000000000000000000000000000..48f638de35557b975d1eb303a16cf42984d1dcdd --- /dev/null +++ b/python/dune/perftool/ufl/rank.py @@ -0,0 +1,15 @@ +from __future__ import absolute_import +from ufl.algorithms import MultiFunction + +class _UFLRank(MultiFunction): + def __call__(self, expr): + return len(MultiFunction.__call__(self, expr)) + + def expr(self, o): + return set(a for op in o.ufl_operands for a in MultiFunction.__call__(self, op)) + + def argument(self, o): + return (o.number(),) + +def ufl_rank(o): + return _UFLRank()(o) \ No newline at end of file diff --git a/python/dune/perftool/ufl/topsum.py b/python/dune/perftool/ufl/topsum.py deleted file mode 100644 index 5f79d6cc88ad4310813a5b78560686ea91497484..0000000000000000000000000000000000000000 --- a/python/dune/perftool/ufl/topsum.py +++ /dev/null @@ -1,25 +0,0 @@ -""" A multifunction for separation of the top sum - -TODO: describe the purpose of topsums. -""" - -from __future__ import absolute_import -from ufl.algorithms import MultiFunction - - -class TopSumSplit(MultiFunction): - """ Split an expression into summands of the top-level sum """ - def __call__(self, o): - self.terms = [] - MultiFunction.__call__(self, o) - return self.terms - - def expr(self, o): - self.terms.append(o) - - def sum(self, o): - for op in o.operands(): - MultiFunction.__call__(self, op) - - def index_sum(self, o): - MultiFunction.__call__(self, o.operands()[0]) diff --git a/python/dune/perftool/ufl/transformations/__init__.py b/python/dune/perftool/ufl/transformations/__init__.py index e723f548563feb4a0c41e128daf3d362d7f2c4b2..9993be27828e1b11e67d8e494fa3ee06d9f56dbd 100644 --- a/python/dune/perftool/ufl/transformations/__init__.py +++ b/python/dune/perftool/ufl/transformations/__init__.py @@ -16,8 +16,7 @@ class TransformationWrapper(object): # Write out a dot file from dune.perftool.options import get_option - # TODO This should be siabled by default! - if get_option("print_transformations", True): + if get_option("print_transformations", False): import os dir = get_option("print_transformations_dir", os.getcwd()) filename = "trafo_{}_{}_{}.dot".format(self.name, str(self.counter).zfill(4), "in" if before else "out") @@ -57,4 +56,18 @@ def ufl_transformation(**kwargs): @ufl_transformation(name="print", printBefore=False) def print_expression(e): - return e \ No newline at end of file + return e + +def transform_integral(integral, trafo): + from ufl import Integral + assert isinstance(integral, Integral) + assert isinstance(trafo, TransformationWrapper) + + return integral.reconstruct(integrand=trafo(integral.integrand())) + +def transform_form(form, trafo): + from ufl import Form + assert isinstance(form, Form) + assert isinstance(trafo, TransformationWrapper) + + return Form([transform_integral(i, trafo) for i in form.integrals()]) diff --git a/python/dune/perftool/ufl/transformations/argument_elimination.py b/python/dune/perftool/ufl/transformations/argument_elimination.py new file mode 100644 index 0000000000000000000000000000000000000000..b2e0afa91d56588c6f016274caebd2b457f26e23 --- /dev/null +++ b/python/dune/perftool/ufl/transformations/argument_elimination.py @@ -0,0 +1,59 @@ +from __future__ import absolute_import +from ufl.algorithms import MultiFunction + + +class EliminateArguments(MultiFunction): + """ Application scope: Take a valid accumulation term and remove the modified argument from it! """ + call = MultiFunction.__call__ + + def __call__(self, o): + from dune.perftool.ufl.modified_terminals import ModifiedArgumentExtractor + from dune.perftool.ufl.rank import ufl_rank + + self.arguments = ModifiedArgumentExtractor()(o) + e = self.call(o) + + # Catch the case that the entire expression vanished! + if e is None: + from ufl.classes import IntValue + e = IntValue(1) + + return (e, self.arguments) + + def expr(self, o): + if o in self.arguments: + return None + else: + # Evaluate the multi function applied to the operands + newop = tuple(self.call(op) for op in o.ufl_operands) + # Find out whether an operand vanished. If so, the class needs special treatment. + if None in newop: + raise NotImplementedError("Operand vanished: {} needs special treatment in EliminateArguments".format(type(o))) + return self.reuse_if_untouched(o, *newop) + + def sum(self, o): + assert len(o.ufl_operands) == 2 + + op0 = self.call(o.ufl_operands[0]) + op1 = self.call(o.ufl_operands[1]) + + if op0 and op1: + return self.reuse_if_untouched(o, op0, op1) + # One term vanished, so there is no sum anymore + else: + if op0 or op1: + # Return the term that did not vanish + return op0 if op0 else op1 + else: + # This entire sum vanished!!! + return None + + # The handler for product is equal to the sum handler + product = sum + + def index_sum(self, o): + op = self.call(o.ufl_operands[0]) + if op: + return self.reuse_if_untouched(o, op, o.ufl_operands[1]) + else: + return None diff --git a/python/dune/perftool/ufl/transformations/extract_accumulation_terms.py b/python/dune/perftool/ufl/transformations/extract_accumulation_terms.py new file mode 100644 index 0000000000000000000000000000000000000000..81feeae68c081e90900e199382d95a3a0c36235f --- /dev/null +++ b/python/dune/perftool/ufl/transformations/extract_accumulation_terms.py @@ -0,0 +1,63 @@ +from __future__ import absolute_import +from ufl.algorithms import MultiFunction + + +class SplitIntoAccumulationTerms(MultiFunction): + """ return a list of tuples of expressions and modified arguments! """ + + call = MultiFunction.__call__ + + def __call__(self, expr): + from dune.perftool.ufl.rank import ufl_rank + self.rank = ufl_rank(expr) + + from dune.perftool.ufl.modified_terminals import ModifiedArgumentExtractor + self.mae = ModifiedArgumentExtractor() + + from dune.perftool.ufl.transformations.argument_elimination import EliminateArguments + self.ea = EliminateArguments() + + # Collect the found terms + self.terms = {} + + # Now, fill the terms dict by applying this multifunction! + self.call(expr) + + from dune.perftool.ufl.flatoperators import construct_binary_operator + from ufl.classes import Sum + return [(construct_binary_operator(t, Sum), a) for a, t in self.terms.items()] + + def expr(self, o): + # This needs to be a valid accumulation term! + assert all(len(self.mae(o, i)) == 1 for i in range(self.rank)) + + expr, arg = self.ea(o) + if arg not in self.terms: + self.terms[arg] = [] + self.terms[arg].append(expr) + + def sum(self, o): + # Check whether this sums contains too many accumulation terms. + if not all(len(self.mae(o, i)) == 1 for i in range(self.rank)): + # This sum is part of a top level sum that separates accumulation terms! + for op in o.ufl_operands: + self.call(op) + else: + # This is a normal sum, we might treat it as any other expression + self.expr(o) + + def index_sum(self, o): + # Check whether this sums contains too many accumulation terms. + if not all(len(self.mae(o, i)) == 1 for i in range(self.rank)): + # This sum is part of a top level sum that separates accumulation terms! + self.call(o.ufl_operands[0]) + else: + # This is a normal sum, we might treat it as any other expression + # TODO we need to eliminate topsum indexsum regardless of the thing being valid + self.expr(o.ufl_operands[0]) + # old code + #self.expr(o) + + +def split_into_accumulation_terms(expr): + return SplitIntoAccumulationTerms()(expr) diff --git a/python/dune/perftool/ufl/transformations/splitarguments.py b/python/dune/perftool/ufl/transformations/splitarguments.py new file mode 100644 index 0000000000000000000000000000000000000000..9304f2ab33f8a56330b574839ac4b2027826ec4f --- /dev/null +++ b/python/dune/perftool/ufl/transformations/splitarguments.py @@ -0,0 +1,46 @@ +from __future__ import absolute_import +from ufl.algorithms import MultiFunction +from dune.perftool.ufl.transformations import ufl_transformation + + +class SplitArguments(MultiFunction): + + call = MultiFunction.__call__ + + def __init__(self): + MultiFunction.__init__(self) + + def __call__(self, o): + from dune.perftool.ufl.modified_terminals import ModifiedArgumentExtractor + self.ae = ModifiedArgumentExtractor() + from dune.perftool.ufl.rank import ufl_rank + self.rank = ufl_rank(o) + # call the actual recursive function + return self.call(o) + + def expr(self, o): + return self.reuse_if_untouched(o, *tuple(self.call(op) for op in o.ufl_operands)) + + def product(self, o): + from dune.perftool.ufl.flatoperators import get_operands, construct_binary_operator + from itertools import product as iterproduct + from ufl.classes import Sum, Product + + # Check whether this product is fine! + if len(self.ae(o)) == self.rank: + return self.reuse_if_untouched(o, *(self.call(op) for op in o.ufl_operands)) + + # It is not, lets switch sums and products! + # First we apply recursively to all + product_operands = [get_operands(self.call(op)) if isinstance(op, Sum) else (self.call(op),) for op in get_operands(o)] + # Multiply all terms by taking the cartesian product of terms + distributive = [f for f in iterproduct(*product_operands)] + # Prepare all sum terms by introducing products + sum_terms = [construct_binary_operator(s, Product) for s in distributive] + # Return the big sum. + return construct_binary_operator(sum_terms, Sum) + + +@ufl_transformation(name="split") +def split_arguments(expr): + return SplitArguments()(expr) diff --git a/python/dune/perftool/ufl/validity.py b/python/dune/perftool/ufl/validity.py index c8c77787de90513bc077d97877b2e212242f8a16..92bee291bad70799ad6b805f614d447d4ca781fe 100644 --- a/python/dune/perftool/ufl/validity.py +++ b/python/dune/perftool/ufl/validity.py @@ -1,26 +1,5 @@ from __future__ import absolute_import - from ufl.algorithms import MultiFunction -from dune.perftool.ufl.modified_terminals import ModifiedTerminalTracker - - -class UFLRank(MultiFunction): - def __call__(self, expr): - return len(MultiFunction.__call__(self, expr)) - - def expr(self, o): - return set(a for op in o.operands() for a in MultiFunction.__call__(self, op)) - - def argument(self, o): - return (o.number(),) - - -class ArgumentExtractor(ModifiedTerminalTracker): - def expr(self, o): - return set(a for op in o.operands() for a in MultiFunction.__call__(self, op)) - - def argument(self, o): - return ((o.number(), self.grad, self.reference_grad, self.restriction),) def check_validity(uflexpr): @@ -33,27 +12,24 @@ def check_validity(uflexpr): """ from ufl import Form from ufl.classes import Expr - from dune.perftool.ufl.topsum import TopSumSplit - - tss = TopSumSplit() - ae = ArgumentExtractor() - def check(term, rank): - if not len(ae(term)) == rank: + def check(term): + try: + from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms + split_into_accumulation_terms(term) + return True + except AssertionError: from dune.perftool.ufl.transformations import print_expression print_expression(term) - raise ValueError("Form not valid for pdelab code generation. Dumped dot file to send to dominic.kempf@iwr.uni-heidelberg.de") + return False if isinstance(uflexpr, Form): - rank = len(uflexpr.arguments()) for integral in uflexpr.integrals(): - for term in tss(integral.integrand()): - check(term, rank) - return + if not check(integral.integrand()): + return False + return True if isinstance(uflexpr, Expr): - rank = UFLRank()(uflexpr) - check(uflexpr, rank) - return + return check(uflexpr) raise TypeError("Unknown object type in check_validity: {}".format(type(uflexpr)))