Merge branch 'feature/refactor-manual-cse' into 'master'

Feature/refactor manual cse

See merge request !114

applications/convection_diffusion/conv_diff_dg.ufl

@@ -16,20 +16,20 @@ n = FacetNormal(cell)('+')
 
 alpha = 3.0
 h_ext = CellVolume(cell) / FacetArea(cell)
-gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+cse_gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
 h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
-gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+cse_gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
 
 theta = -1.0
 
 r = (inner(A*grad(u), grad(v)) + (c*u-f)*v)*dx \
   + inner(n, A*avg(grad(u)))*jump(v)*dS \
-  + gamma_int*jump(u)*jump(v)*dS \
+  + cse_gamma_int*jump(u)*jump(v)*dS \
   - theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
   - inner(n, A*grad(u))*v*ds \
-  + gamma_ext*u*v*ds \
+  + cse_gamma_ext*u*v*ds \
   + theta*u*inner(A*grad(v), n)*ds \
   - theta*g*inner(A*grad(v), n)*ds \
-  - gamma_ext*g*v*ds
+  - cse_gamma_ext*g*v*ds
 
 forms = [r]

applications/poisson_dg/poisson_dg.ufl

@@ -12,21 +12,21 @@ n = FacetNormal(cell)('+')
 
 alpha = 1.0
 h_ext = CellVolume(cell) / FacetArea(cell)
-gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+cse_gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
 h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
-gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+cse_gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
 
 theta = -1.0
 
 r = inner(grad(u), grad(v))*dx \
   + inner(n, avg(grad(u)))*jump(v)*dS \
-  + gamma_int*jump(u)*jump(v)*dS \
+  + cse_gamma_int*jump(u)*jump(v)*dS \
   - theta*jump(u)*inner(avg(grad(v)), n)*dS \
   - inner(n, grad(u))*v*ds \
-  + gamma_ext*u*v*ds \
+  + cse_gamma_ext*u*v*ds \
   + theta*u*inner(grad(v), n)*ds \
   - f*v*dx \
   - theta*g*inner(grad(v), n)*ds \
-  - gamma_ext*g*v*ds
+  - cse_gamma_ext*g*v*ds
 
 forms = [r]

python/dune/perftool/compile.py

@@ -12,7 +12,6 @@ from ufl.algorithms.formfiles import interpret_ufl_namespace
 
 from dune.perftool.generation import (delete_cache_items,
                                       global_context,
-                                      subst_rule,
                                       )
 from dune.perftool.interactive import start_interactive_session
 from dune.perftool.options import get_option

python/dune/perftool/generation/__init__.py

@@ -33,6 +33,7 @@ from dune.perftool.generation.cpp import (base_class,
 from dune.perftool.generation.loopy import (barrier,
                                             built_instruction,
                                             constantarg,
+                                            construct_subst_rule,
                                             domain,
                                             function_mangler,
                                             get_temporary_name,

python/dune/perftool/generation/loopy.py

@@ -201,15 +201,29 @@ def loopy_class_member(name, classtag=None, potentially_vectorized=False, **kwar
     return name
 
 
-@generator_factory(item_tags=("substrule",),
+@generator_factory(item_tags=("substrule_name",),
+                   context_tags="kernel",
+                   cache_key_generator=lambda e, n: e)
+def _substrule_name(expr, name):
+    return name
+
+
+@generator_factory(item_tags=("substrule_impl",),
                    context_tags="kernel",
-                   cache_key_generator=lambda e, r: e,
+                   cache_key_generator=lambda n, e, **ex: e,
                    )
-def subst_rule(expr, rule):
-    return rule
+def subst_rule(name, expr, exists=False):
+    _substrule_name(expr, name)
+    return exists
+
 
+def set_subst_rule(name, expr):
+    subst_rule(name, expr, exists=True)
 
-def set_subst_rule(name, expr, visitor):
-    rule = lp.SubstitutionRule(name, (), visitor(expr))
-    subst_rule._memoize_cache = {k: v for k, v in subst_rule._memoize_cache.items() if v is not None}
-    return subst_rule(expr, rule)
+
+@generator_factory(item_tags=("substrule",),
+                   context_tags="kernel")
+def construct_subst_rule(expr, visitor):
+    name = _substrule_name(expr, None)
+    assert name
+    return lp.SubstitutionRule(name, (), visitor(expr))

python/dune/perftool/loopy/transformations/collect_rotate.py

@@ -194,7 +194,7 @@ def collect_vector_data_rotate(knl):
                     replacemap_vec[expr] = prim.Subscript(prim.Variable(get_vector_view_name(quantity)),
                                                           (vector_indices.get(vec_size) + last_index, prim.Variable(new_iname)),
                                                           )
-            else:
+            elif all(get_pymbolic_tag(expr) == 'sumfac' for expr in quantities[quantity]):
                 # Add a vector view to this quantity
                 expr, = quantities[quantity]
                 knl = add_vector_view(knl, quantity, flatview=True)

python/dune/perftool/pdelab/localoperator.py

@@ -419,16 +419,16 @@ def visit_integrals(integrals):
             data = get_global_context_value("data")
             for name, expr in data.object_by_name.items():
                 if name.startswith("cse"):
-                    set_subst_rule(name, expr, visitor)
+                    set_subst_rule(name, expr)
 
             # Ensure CSE on detjac * quadrature weight
             domain = term.argument.argexpr.ufl_domain()
             if measure == "cell":
-                set_subst_rule("integration_factor_cell1", uc.QuadratureWeight(domain)*uc.Abs(uc.JacobianDeterminant(domain)), visitor)
-                set_subst_rule("integration_factor_cell2", uc.Abs(uc.JacobianDeterminant(domain))*uc.QuadratureWeight(domain), visitor)
+                set_subst_rule("integration_factor_cell1", uc.QuadratureWeight(domain)*uc.Abs(uc.JacobianDeterminant(domain)))
+                set_subst_rule("integration_factor_cell2", uc.Abs(uc.JacobianDeterminant(domain))*uc.QuadratureWeight(domain))
             else:
-                set_subst_rule("integration_factor_facet1", uc.FacetJacobianDeterminant(domain)*uc.QuadratureWeight(domain), visitor)
-                set_subst_rule("integration_factor_facet2", uc.QuadratureWeight(domain)*uc.FacetJacobianDeterminant(domain), visitor)
+                set_subst_rule("integration_factor_facet1", uc.FacetJacobianDeterminant(domain)*uc.QuadratureWeight(domain))
+                set_subst_rule("integration_factor_facet2", uc.QuadratureWeight(domain)*uc.FacetJacobianDeterminant(domain))
 
             get_backend(interface="accum_insn")(visitor, term, measure, subdomain_id)
 

python/dune/perftool/sumfact/amatrix.py

@@ -54,7 +54,7 @@ class AMatrix(ImmutableRecord):
         return False
 
     def output_to_pymbolic(self, name):
-        return prim.Variable(name)
+        return lp.TaggedVariable(name, "sumfac")
 
 
 class LargeAMatrix(ImmutableRecord):

python/dune/perftool/ufl/visitor.py

@@ -3,7 +3,8 @@ This module defines the main visitor algorithm transforming ufl expressions
 to pymbolic and loopy.
 """
 from dune.perftool.error import PerftoolUFLError
-from dune.perftool.generation import (domain,
+from dune.perftool.generation import (construct_subst_rule,
+                                      domain,
                                       get_global_context_value,
                                       subst_rule,
                                       )
@@ -38,10 +39,11 @@ import pymbolic.primitives as prim
 
 
 class UFL2LoopyVisitor(ModifiedTerminalTracker):
-    def __init__(self, interface, measure, dimension_indices):
+    def __init__(self, interface, measure, dimension_indices, donot_check_substrules=None):
         self.interface = interface
         self.measure = measure
         self.dimension_indices = dimension_indices
+        self.donot_check_substrules = donot_check_substrules
 
         # Call base class constructors
         super(UFL2LoopyVisitor, self).__init__()
@@ -55,8 +57,14 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         return self.call(o)
 
     def call(self, o):
-        rule = subst_rule(o, None)
-        if rule:
+        if o != self.donot_check_substrules and subst_rule(None, o):
+            rule = construct_subst_rule(o,
+                                        type(self)(self.interface,
+                                                   self.measure,
+                                                   self.dimension_indices,
+                                                   donot_check_substrules=o,
+                                                   )
+                                        )
             return prim.Call(prim.Variable(rule.name), ())
         else:
             return MultiFunction.__call__(self, o)