Merge branch 'feature/sumfact-stokes-dg-numdiff' into 'master'

[Feature] Stokes DG (numdiff) with sum factorization

See merge request !146

python/dune/perftool/pdelab/localoperator.py

@@ -281,7 +281,7 @@ def boundary_predicates(expr, measure, subdomain_id, visitor):
 
         from ufl.classes import Expr
         if isinstance(subdomain_data, Expr):
-            cond = visitor(subdomain_data)
+            cond = visitor(subdomain_data, do_predicates=True)
         else:
             # Determine the name of the parameter function
             cond = get_global_context_value("data").object_names[id(subdomain_data)]

python/dune/perftool/sumfact/__init__.py

@@ -61,6 +61,6 @@ class SumFactInterface(PDELabInterface):
 
     def pymbolic_spatial_coordinate(self):
         import dune.perftool.sumfact.geometry
-        ret, indices = get_backend(interface="spatial_coordinate", selector=option_switch("diagonal_transformation_matrix"))(self.visitor.indices)
+        ret, indices = get_backend(interface="spatial_coordinate", selector=option_switch("diagonal_transformation_matrix"))(self.visitor.indices, self.visitor.do_predicates)
         self.visitor.indices = indices
         return ret

python/dune/perftool/sumfact/accumulation.py

@@ -68,8 +68,42 @@ class AlreadyAssembledInput(SumfactKernelInputBase):
     def __eq__(self, other):
         return type(self) == type(other) and self.index == other.index
 
+    def __repr__(self):
+        return "AlreadyAssembledInput({})".format(self.index)
+
     def __hash__(self):
-        return 0
+        return hash(self.index)
+
+
+def _dimension_gradient_indices(argument):
+    """Return dimension and gradient index of test function argument
+
+    Dimension index corresponds to the dimensio for vector valued
+    functions (eg v in Stokes). The gradient index is the direction of
+    the derivative (eg in \Delat u in poisson or \grad u in Stoks).
+    """
+    grad_index = None
+    dimension_index = None
+
+    # If this is a gradient the last index is the grad_index
+    if argument.reference_grad:
+        grad_index = argument.expr.ufl_operands[1][-1]._value
+
+    # If this argument has indices there could be a dimension_index
+    if isinstance(argument.expr, uc.Indexed):
+        # More than two indices not supported
+        if len(argument.expr.ufl_operands[1]) > 2:
+            assert False
+        # For two indices the first is the dimension index
+        if len(argument.expr.ufl_operands[1]) == 2:
+            assert grad_index is not None
+            dimension_index = argument.expr.ufl_operands[1].indices()[0]._value
+        # For there is no gradient index we should have only one index, the dimension index
+        if not argument.reference_grad:
+            assert len(argument.expr.ufl_operands[1]) == 1
+            dimension_index = argument.expr.ufl_operands[1].indices()[0]._value
+
+    return dimension_index, grad_index
 
 
 @backend(interface="accum_insn", name="sumfact")
@@ -77,37 +111,35 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
     # When doing sum factorization we want to split the test function
     assert(accterm.argument.expr is not None)
 
+    # Get dimension and gradient index
+    dimension_index, grad_index = _dimension_gradient_indices(accterm.argument)
+    accum_index = (dimension_index,)
+
     # Do the tree traversal to get a pymbolic expression representing this expression
     pymbolic_expr = visitor(accterm.term)
     if pymbolic_expr == 0:
         return
 
+    # Number of basis functions
     dim = world_dimension()
     mod_arg_expr = accterm.argument.expr
-    from ufl.classes import FunctionView, Argument
-    while (not isinstance(mod_arg_expr, FunctionView)) and (not isinstance(mod_arg_expr, Argument)):
+    while (not isinstance(mod_arg_expr, uc.FunctionView)) and (not isinstance(mod_arg_expr, uc.Argument)):
         mod_arg_expr = mod_arg_expr.ufl_operands[0]
+
     degree = mod_arg_expr.ufl_element()._degree
     basis_size = degree + 1
 
-    # Extract index information
-    grad_index = None
-    if accterm.argument.reference_grad:
-        grad_index = accterm.argument.expr.ufl_operands[1][0]._value
-
-    accum_index = None
-    if visitor.indices:
-        accum_index = visitor.indices[0]
-    if accterm.argument.index and len(accterm.argument.index) > 1:
-        accum_index = accterm.argument.index[1]._value
-
     jacobian_inames = tuple()
     if accterm.is_jacobian:
         jacobian_inames = visitor.inames
 
+    # Only accumulate boundary conditions on parts where it is defined
+    from dune.perftool.pdelab.localoperator import boundary_predicates
+    predicates = boundary_predicates(accterm.term, measure, subdomain_id, visitor)
+
     def emit_sumfact_kernel(restriction, insn_dep):
-        test_lfs = determine_accumulation_space(accterm.argument.expr, 0, measure, idims=(accum_index,))
-        ansatz_lfs = determine_accumulation_space(accterm.term, 1, measure, idims=(accum_index,))
+        test_lfs = determine_accumulation_space(accterm.argument.expr, 0, measure, idims=accum_index)
+        ansatz_lfs = determine_accumulation_space(accterm.term, 1, measure, idims=accum_index)
         accum = name_accumulation_variable(test_lfs.get_restriction() + ansatz_lfs.get_restriction())
 
         # Construct the matrix sequence for this sum factorization
@@ -126,6 +158,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
                            accumvar=accum,
                            within_inames=jacobian_inames,
                            input=AlreadyAssembledInput(index=accum_index),
+                           coeff_func_index=dimension_index,
                            )
 
         from dune.perftool.sumfact.vectorization import attach_vectorization_info
@@ -152,7 +185,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
                         expression=0,
                         forced_iname_deps=frozenset(quadrature_inames() + jacobian_inames),
                         forced_iname_deps_is_final=True,
-                        tags=frozenset(["quadvec", "gradvec"])
+                        tags=frozenset(["quadvec", "gradvec"]),
                         )
 
         # Write timing stuff for jacobian (for alpha methods it is done at the end of stage 1)
@@ -180,7 +213,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
                                   forced_iname_deps=frozenset(quadrature_inames() + jacobian_inames),
                                   forced_iname_deps_is_final=True,
                                   tags=frozenset({"quadvec"}).union(vectag),
-                                  depends_on=frozenset({deps}).union(timer_dep).union(frozenset({lp.match.Tagged("sumfact_stage1")}))
+                                  depends_on=frozenset({deps}).union(timer_dep).union(frozenset({lp.match.Tagged("sumfact_stage1")})),
                                   )
 
         if insn_dep is None:
@@ -241,6 +274,7 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
                         forced_iname_deps=frozenset(inames + vecinames + jacobian_inames),
                         forced_iname_deps_is_final=True,
                         depends_on=insn_dep,
+                        predicates=predicates
                         )
 
         # Mark the transformation that moves the quadrature loop

python/dune/perftool/sumfact/basis.py

@@ -150,7 +150,7 @@ def pymbolic_coefficient_gradient(element, restriction, component, coeff_func, v
         coeff_func_index = None
         grad_index, = visitor_indices
     else:
-        grad_index, coeff_func_index = visitor_indices
+        coeff_func_index, grad_index = visitor_indices
 
     # Construct the matrix sequence for this sum factorization
     matrix_sequence = construct_basis_matrix_sequence(derivative=grad_index,
@@ -183,7 +183,6 @@ def pymbolic_coefficient_gradient(element, restriction, component, coeff_func, v
 
 @kernel_cached
 def pymbolic_coefficient(element, restriction, component, coeff_func, visitor_indices):
-    assert visitor_indices is None
     # Basis functions per direction
     basis_size = _basis_functions_per_direction(element, component)
 
@@ -192,7 +191,13 @@ def pymbolic_coefficient(element, restriction, component, coeff_func, visitor_in
                                                       facemod=get_facemod(restriction),
                                                       basis_size=basis_size)
 
+    coeff_func_index = None
+    if visitor_indices:
+        assert len(visitor_indices) == 1
+        coeff_func_index = visitor_indices[0]
+
     inp = LFSSumfactKernelInput(coeff_func=coeff_func,
+                                coeff_func_index=coeff_func_index,
                                 element=element,
                                 component=component,
                                 restriction=restriction,

python/dune/perftool/sumfact/geometry.py

@@ -68,7 +68,7 @@ class GeoCornersInput(SumfactKernelInputBase, ImmutableRecord):
 
 @kernel_cached
 @backend(interface="spatial_coordinate", name="default")
-def pymbolic_spatial_coordinate_multilinear(visitor_indices):
+def pymbolic_spatial_coordinate_multilinear(visitor_indices, do_predicates):
     assert len(visitor_indices) == 1
 
     # Construct the matrix sequence for the evaluation of the global coordinate.
@@ -127,7 +127,7 @@ def name_meshwidth():
 
 @kernel_cached
 @backend(interface="spatial_coordinate", name="diagonal_transformation_matrix")
-def pymbolic_spatial_coordinate_axiparallel(visitor_indices):
+def pymbolic_spatial_coordinate_axiparallel(visitor_indices, do_predicates):
     assert len(visitor_indices) == 1
     index, = visitor_indices
 
@@ -143,7 +143,13 @@ def pymbolic_spatial_coordinate_axiparallel(visitor_indices):
     lowcorner = name_lowerleft_corner()
     meshwidth = name_meshwidth()
 
-    if index == face:
+    # If we have to decide which boundary condition to take for this
+    # intersection we always take the boundary condition of the center
+    # of the intersection. We assume that there are no cells with more
+    # than one boundary condition.
+    if do_predicates:
+        x = 0.5
+    elif index == face:
         x = 0
     else:
         iindex = index

python/dune/perftool/sumfact/symbolic.py

@@ -38,6 +38,7 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
                  insn_dep=frozenset(),
                  input=None,
                  accumvar=None,
+                 coeff_func_index=None,
                  ):
         """Create a sum factorization kernel
 
@@ -91,23 +92,18 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
             pre-initialized with the input or you have to provide
             direct_input (FastDGGridOperator).
         stage: 1 or 3
-        insn_dep: An instruction ID that the first issued instruction
-            should depend upon. All following ones will depend on each
-            other.
-        within_inames: Instructions will be executed within those
-            inames (needed for stage 3 in jacobians).
         preferred_position: Will be used in the dry run to order kernels
             when doing vectorization e.g. (dx u,dy u,dz u, u).
         restriction: Restriction for faces values.
-        padding: a set of slots in the input temporary to be padded
-        index: The slot in the input temporary dedicated to this kernel
-        coeff_func: The function to call to get the input coefficient
+        within_inames: Instructions will be executed within those
+            inames (needed for stage 3 in jacobians).
+        insn_dep: An instruction ID that the first issued instruction
+            should depend upon. All following ones will depend on each
+            other.
+        input: An SumfactKernelInputBase instance describing the input of the kernel
+        accumvar: The accumulation variable to accumulate into
         coeff_func_index: Index to get the right input coefficient
             (needed for systems of PDEs)
-        element: The UFL element
-        component: The treepath to the correct component of above element
-        accumvar: The accumulation variable to accumulate into
-        input: An SumfactKernelInputBase instance describing the input of the kernel
         """
         # Assert the inputs!
         assert isinstance(matrix_sequence, tuple)
@@ -162,7 +158,7 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
         Any two sum factorization kernels having the same cache_key
         are realized simulatenously!
         """
-        return (self.matrix_sequence, self.restriction, self.stage, self.buffer)
+        return (self.matrix_sequence, self.restriction, self.stage, self.buffer, self.coeff_func_index)
 
     @property
     def input_key(self):

python/dune/perftool/ufl/visitor.py

@@ -44,12 +44,12 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
         # Call base class constructors
         super(UFL2LoopyVisitor, self).__init__()
 
-    def __call__(self, o):
+    def __call__(self, o, do_predicates=False):
         # Reset some state variables that are reinitialized for each accumulation term
         self.indices = None
         self._indices_backup = []
         self.inames = ()
-
+        self.do_predicates = do_predicates
         return self.call(o)
 
     call = MultiFunction.__call__

test/stokes/stokes_dg_quadrilateral.ufl

@@ -20,18 +20,18 @@ eps = -1.0
 sigma = 1.0
 
 r = inner(grad(u), grad(v))*dx \
+  - p*div(v)*dx \
+  - q*div(u)*dx \
   + inner(avg(grad(u))*n, jump(v))*dS \
-  - eps * inner(avg(grad(v))*n, jump(u))*dS \
-  - inner(grad(u)*n, v)*ds \
-  + eps * inner(grad(v)*n, u-g_v)*ds \
   + sigma * inner(jump(u), jump(v))*dS \
-  + sigma * inner(u-g_v, v)*ds \
-  - p*div(v)*dx \
+  - eps * inner(avg(grad(v))*n, jump(u))*dS \
   - avg(p)*inner(jump(v), n)*dS \
-  + p*inner(v, n)*ds \
-  - q*div(u)*dx \
   - avg(q)*inner(jump(u), n)*dS \
+  - inner(grad(u)*n, v)*ds \
+  + p*inner(v, n)*ds \
   + q*inner(u, n)*ds \
+  + eps * inner(grad(v)*n, u-g_v)*ds \
+  + sigma * inner(u-g_v, v)*ds \
   - q*inner(g_v, n)*ds
 
 forms = [r]

test/sumfact/stokes/CMakeLists.txt

@@ -2,3 +2,8 @@ dune_add_formcompiler_system_test(UFLFILE stokes.ufl
                                   BASENAME sumfact_stokes
                                   INIFILE stokes.mini
                                   )
+
+dune_add_formcompiler_system_test(UFLFILE stokes_dg.ufl
+                                  BASENAME sumfact_stokes_dg
+                                  INIFILE stokes_dg.mini
+                                  )

test/sumfact/stokes/stokes_dg.mini

+__name = sumfact_stokes_dg_{__exec_suffix}
+
+__exec_suffix = symdiff, numdiff | expand num
+
+cells = 8 8
+extension = 1. 1.
+printmatrix = false
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-8
+sumfact = 1
+
+# Disable symdiff test for now
+{__exec_suffix} == symdiff | exclude

test/sumfact/stokes/stokes_dg.ufl

+cell = quadrilateral
+
+x = SpatialCoordinate(cell)
+g_v = as_vector((4*x[1]*(1.-x[1]), 0.0))
+g_p = 8*(1.-x[0])
+g = (g_v, g_p)
+bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
+
+P2 = VectorElement("DG", cell, 2)
+P1 = FiniteElement("DG", cell, 1)
+TH = P2 * P1
+
+v, q = TestFunctions(TH)
+u, p = TrialFunctions(TH)
+
+ds = ds(subdomain_id=1, subdomain_data=bctype)
+
+n = FacetNormal(cell)('+')
+eps = -1.0
+sigma = 1.0
+
+r = inner(grad(u), grad(v))*dx \
+  - p*div(v)*dx \
+  - q*div(u)*dx \
+  + inner(avg(grad(u))*n, jump(v))*dS \
+  + sigma * inner(jump(u), jump(v))*dS \
+  - eps * inner(avg(grad(v))*n, jump(u))*dS \
+  - avg(p)*inner(jump(v), n)*dS \
+  - avg(q)*inner(jump(u), n)*dS \
+  - inner(grad(u)*n, v)*ds \
+  + p*inner(v, n)*ds \
+  + q*inner(u, n)*ds \
+  + eps * inner(grad(v)*n, u-g_v)*ds \
+  + sigma * inner(u-g_v, v)*ds \
+  - q*inner(g_v, n)*ds
+
+forms = [r]