WIP: Working on sumfactorization for unstructered grids

Not yet implemented:

- Vectorization of quadrature loop
- Anisotropic amount of quadrature points per direction

python/dune/perftool/pdelab/geometry.py

@@ -253,11 +253,11 @@ def local_dimension():
 def evaluate_unit_outer_normal(name):
     ig = name_intersection_geometry_wrapper()
     qp = get_backend("quad_pos")()
-    return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
-                               assignees=frozenset({name}),
-                               read_variables=frozenset({get_pymbolic_basename(qp)}),
-                               depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
-                               )
+    return get_backend("quadrature_preamble")("{} = {}.unitOuterNormal({});".format(name, ig, qp),
+                                              assignees=frozenset({name}),
+                                              read_variables=frozenset({get_pymbolic_basename(qp)}),
+                                              depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
+                                              )
 
 
 @preamble

python/dune/perftool/sumfact/__init__.py

@@ -78,8 +78,10 @@ class SumFactInterface(PDELabInterface):
         return ret
 
     def pymbolic_unit_outer_normal(self):
-        from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
-        ret, indices = pymbolic_unit_outer_normal(self.visitor.indices)
+        from dune.perftool.generation import global_context
+        with global_context(visitor=self.visitor):
+            from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
+            ret, indices = pymbolic_unit_outer_normal(self.visitor.indices)
         self.visitor.indices = indices
         return ret
 
@@ -90,9 +92,23 @@ class SumFactInterface(PDELabInterface):
         return ret
 
     def pymbolic_facet_jacobian_determinant(self):
-        from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
-        return pymbolic_facet_jacobian_determinant()
+        from dune.perftool.generation import global_context
+        with global_context(visitor=self.visitor):
+            from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
+            return pymbolic_facet_jacobian_determinant()
 
     def pymbolic_facet_area(self):
         from dune.perftool.sumfact.geometry import pymbolic_facet_area
         return pymbolic_facet_area()
+
+    def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
+        from dune.perftool.pdelab.geometry import pymbolic_jacobian_inverse_transposed
+        from dune.perftool.generation import global_context
+        with global_context(visitor=self.visitor):
+            return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+
+    def pymbolic_jacobian_determinant(self):
+        from dune.perftool.pdelab.geometry import pymbolic_jacobian_determinant
+        from dune.perftool.generation import global_context
+        with global_context(visitor=self.visitor):
+            return pymbolic_jacobian_determinant()

python/dune/perftool/sumfact/geometry.py

@@ -39,7 +39,17 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
     def __init__(self, dir):
         ImmutableRecord.__init__(self, dir=dir)
 
-    def realize(self, sf, index, insn_dep):
+    @property
+    def stage(self):
+        return 1
+
+    def __repr__(self):
+        return ImmutableRecord.__repr__(self)
+
+    def realize(self, sf, insn_dep):
+        # TODO: When we use vectorization this should be the offset
+        index = 0
+
         from dune.perftool.sumfact.realization import name_buffer_storage
         name = "input_{}".format(sf.buffer)
         temporary_variable(name,
@@ -64,11 +74,13 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
                                                           self.dir,
                                                           )
 
-        instruction(code=code,
-                    within_inames=frozenset({ciname}),
-                    assignees=(name,),
-                    tags=frozenset({"sumfact_stage{}".format(sf.stage)}),
-                    )
+        insn = instruction(code=code,
+                           within_inames=frozenset({ciname}),
+                           assignees=(name,),
+                           tags=frozenset({"sumfact_stage{}".format(sf.stage)}),
+                           )
+
+        return insn_dep.union(frozenset({insn}))
 
 
 @backend(interface="spatial_coordinate", name="default")

python/dune/perftool/sumfact/quadrature.py

@@ -2,13 +2,15 @@ from dune.perftool.generation import (backend,
                                       domain,
                                       function_mangler,
                                       get_global_context_value,
+                                      globalarg,
                                       iname,
                                       instruction,
                                       kernel_cached,
                                       loopy_class_member,
+                                      preamble,
                                       temporary_variable,
                                       )
-
+from dune.perftool.sumfact.switch import get_facedir
 from dune.perftool.sumfact.tabulation import (quadrature_points_per_direction,
                                               local_quadrature_points_per_direction,
                                               name_oned_quadrature_points,
@@ -78,7 +80,6 @@ def pymbolic_base_weight():
     return 1.0
 
 
-@backend(interface="quad_inames", name="sumfact")
 @kernel_cached
 def quadrature_inames(element):
     if element is None:
@@ -179,9 +180,10 @@ def quadrature_weight(visitor):
     return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
 
 
-def define_quadrature_position(name, index):
+def define_quadrature_position(name, local_index):
     qps_per_dir = quadrature_points_per_direction()
-    qp_bound = qps_per_dir[index]
+    local_qps_per_dir = local_quadrature_points_per_direction()
+    qp_bound = local_qps_per_dir[index]
     expression = Subscript(Variable(name_oned_quadrature_points(qp_bound)),
                            (Variable(quadrature_inames()[index]),))
     instruction(expression=expression,
@@ -192,20 +194,29 @@ def define_quadrature_position(name, index):
                 )
 
 
-def pymbolic_indexed_quadrature_position(index, visitor):
+def pymbolic_indexed_quadrature_position(local_index, visitor):
+    """Return the value of the indexed local quadrature points
+
+    Arguments:
+    ----------
+
+    local_index: Local index to the quadrature point. Local means that on the
+        face of a 3D Element the index always will be 0 or 1 but never 2.
+    visitor: UFL tree visitor
+    """
     # Return the non-precomputed version
     if not get_form_option("precompute_quadrature_info"):
         name = 'pos'
         temporary_variable(name, shape=(local_dimension(),), shape_impl=("fv",))
-        define_quadrature_position(name, index)
-        return prim.Subscript(prim.Variable(name), (index,))
+        define_quadrature_position(name, local_index)
+        return prim.Subscript(prim.Variable(name), (local_index,))
 
-    assert isinstance(index, int)
+    assert isinstance(local_index, int)
     dim = local_dimension()
     qps_per_dir = quadrature_points_per_direction()
     local_qps_per_dir = local_quadrature_points_per_direction()
     local_qps_per_dir_str = '_'.join(map(str, local_qps_per_dir))
-    name = "quad_points_dim{}_num{}_dir{}".format(dim, local_qps_per_dir_str, index)
+    name = "quad_points_dim{}_num{}_localdir{}".format(dim, local_qps_per_dir_str, local_index)
 
     loopy_class_member(name,
                        shape=local_qps_per_dir,
@@ -218,8 +229,8 @@ def pymbolic_indexed_quadrature_position(index, visitor):
     # Precompute it in the constructor
     inames = constructor_quadrature_inames(dim)
     assignee = prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames))
-    expression = Subscript(Variable(name_oned_quadrature_points(local_qps_per_dir[index])),
-                           (prim.Variable(inames[index])))
+    expression = Subscript(Variable(name_oned_quadrature_points(local_qps_per_dir[local_index])),
+                           (prim.Variable(inames[local_index])))
     instruction(assignee=assignee,
                 expression=expression,
                 within_inames=frozenset(inames),
@@ -234,10 +245,79 @@ def pymbolic_indexed_quadrature_position(index, visitor):
     return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
 
 
+def _additional_inames(visitor):
+    info = visitor.current_info[1]
+    if info is None:
+        element = None
+    else:
+        element = info.element
+
+    if element is not None:
+        from dune.perftool.pdelab.restriction import Restriction
+        restriction = Restriction.NONE
+        if get_global_context_value("integral_type") == "interior_facet":
+            restriction = Restriction.POSITIVE
+        from dune.perftool.sumfact.basis import lfs_inames
+        lfs_inames = lfs_inames(element, restriction)
+        return lfs_inames
+    else:
+        return ()
+
+
+@backend(interface="quadrature_preamble", name="sumfact")
+def quadrature_preamble(code, **kw):
+    visitor = get_global_context_value('visitor')
+
+    info = visitor.current_info[1]
+    if info is None:
+        element = None
+    else:
+        element = info.element
+
+    inames = quadrature_inames(element)
+    inames = inames + _additional_inames(visitor)
+
+    return instruction(inames=inames, code=code, **kw)
+
+
+@preamble
+def define_quadrature_point(name):
+    local_dim = local_dimension()
+    return "Dune::FieldVector<double, {}> {};".format(local_dim, name)
+
+
+def name_quadrature_point():
+    name = "quadrature_point"
+    define_quadrature_point(name)
+    return name
+
+
+@backend(interface="quad_pos", name="sumfact")
+def pymbolic_quadrature_position():
+    visitor = get_global_context_value('visitor')
+
+    try:
+        info = visitor.current_info[1]
+    except:
+        from pudb import set_trace; set_trace()
+    if info is None:
+        element = None
+    else:
+        element = info.element
+
+    quad_point = name_quadrature_point()
+    inames = quadrature_inames(element) + _additional_inames(visitor)
+    globalarg(quad_point, shape=(local_dimension(),))
+    facedir = get_facedir(visitor.restriction)
+    for i in range(local_dimension()):
+        expression = pymbolic_indexed_quadrature_position(i, visitor)
+        instruction(assignee=prim.Subscript(prim.Variable(quad_point), i),
+                    expression=expression,
+                    within_inames=frozenset(inames))
+    return Variable(quad_point)
+
+
 @backend(interface="qp_in_cell", name="sumfact")
 def pymbolic_quadrature_position_in_cell(restriction):
-    # TODO: This code path is broken at the moment.
-    assert False
-
-    from dune.perftool.pdelab.geometry import to_cell_coordinates
-    return to_cell_coordinates(pymbolic_quadrature_position(), restriction)
+    from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
+    return pymbolic_quadrature_position_in_cell(restriction)

python/dune/perftool/sumfact/switch.py

@@ -10,7 +10,7 @@ from dune.perftool.pdelab.signatures import (assembly_routine_args,
                                              assembly_routine_signature,
                                              kernel_name,
                                              )
-from dune.perftool.options import get_form_option
+from dune.perftool.options import get_form_option, get_option
 from dune.perftool.cgen.clazz import ClassMember
 
 

python/dune/perftool/sumfact/tabulation.py

@@ -20,6 +20,7 @@ from dune.perftool.generation import (class_member,
                                       )
 from dune.perftool.loopy.target import dtype_floatingpoint
 from dune.perftool.loopy.vcl import ExplicitVCLCast, get_vcl_type_size
+from dune.perftool.options import get_option
 from dune.perftool.pdelab.localoperator import (name_domain_field,
                                                 lop_template_range_field,
                                                 )
@@ -273,8 +274,13 @@ def local_quadrature_points_per_direction():
     qps_per_dir = quadrature_points_per_direction()
     if local_dimension() != world_dimension():
         facedir = get_global_context_value('facedir_s')
-        assert(get_global_context_value('facedir_s') == get_global_context_value('facedir_n') or
-               get_global_context_value('integral_type') == 'exterior_facet')
+        if not get_option("grid_unstructured"):
+            assert(get_global_context_value('facedir_s') == get_global_context_value('facedir_n') or
+                   get_global_context_value('integral_type') == 'exterior_facet')
+        if get_option("grid_unstructured"):
+            # For unstructured grids the amount of quadrature points could be different for
+            # self and neighbor. For now assert that this case is not happining.
+            assert len(set(qps_per_dir))==1
         qps_per_dir = qps_per_dir[:facedir] + qps_per_dir[facedir + 1:]
     return qps_per_dir
 

test/sumfact/poisson/poisson_2d_unstructured.mini

@@ -6,6 +6,11 @@ quadvec_suffix = quadvec, nonquadvec | expand quad
 gradvec_suffix = gradvec, nongradvec | expand grad
 deg_suffix = deg{formcompiler.ufl_variants.degree}
 
+{deg_suffix} == deg1 | exclude
+{diff_suffix} == numdiff | exclude
+{quadvec_suffix} == quadvec | exclude
+
+
 lowerleft = 0.0 0.0
 upperright = 1.0 1.0
 elements = 16 16
@@ -25,9 +30,5 @@ sumfact = 1
 vectorization_quadloop = 1, 0 | expand quad
 vectorization_strategy = explicit, none | expand grad
 
-# palpo TODO
-# diagonal_transformation_matrix = 1
-# constant_transformation_matrix = 1
-
 [formcompiler.ufl_variants]
 degree = 1, 2 | expand deg

test/sumfact/poisson/poisson_dg_2d_unstructured.mini

@@ -6,6 +6,11 @@ quadvec_suffix = quadvec, nonquadvec | expand quad
 gradvec_suffix = gradvec, nongradvec | expand grad
 deg_suffix = deg{formcompiler.ufl_variants.degree}
 
+{deg_suffix} == deg1 | exclude
+{diff_suffix} == numdiff | exclude
+{quadvec_suffix} == quadvec | exclude
+
+
 lowerleft = 0.0 0.0
 upperright = 1.0 1.0
 elements = 16 16
@@ -24,9 +29,6 @@ numerical_jacobian = 1, 0 | expand num
 sumfact = 1
 vectorization_quadloop = 1, 0 | expand quad
 vectorization_strategy = explicit, none | expand grad
-# palpo TODO
-diagonal_transformation_matrix = 1
-constant_transformation_matrix = 1
 
 [formcompiler.ufl_variants]
 degree = 1, 2 | expand deg