Also precompute quadrature positions

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

@@ -244,7 +244,7 @@ def collect_vector_data_rotate(knl):
                 shape=(ceildiv(product(s for s in arg.shape), vec_size), vec_size)
                 name = loopy_class_member(quantity,
                                           shape=shape,
-                                          dim_tags="c,vec",
+                                          dim_tags="f,vec",
                                           potentially_vectorized=True,
                                           classtag="operator",
                                           dtype=np.float64,

python/dune/perftool/sumfact/geometry.py

@@ -149,6 +149,6 @@ def pymbolic_spatial_coordinate_axiparallel(visitor_indices):
         if face is not None and index > face:
             iindex = iindex - 1
         from dune.perftool.sumfact.quadrature import pymbolic_quadrature_position
-        x = prim.Subscript(pymbolic_quadrature_position(), (iindex,))
+        x = pymbolic_quadrature_position(iindex)
 
     return prim.Subscript(prim.Variable(lowcorner), (index,)) + x * prim.Subscript(prim.Variable(meshwidth), (index,)), None

python/dune/perftool/sumfact/quadrature.py

@@ -93,7 +93,8 @@ def constructor_quad_iname(name, d, bound):
     return name
 
 
-def constructor_quadrature_inames(name):
+def constructor_quadrature_inames(dim, num1d):
+    name = "quadiname_dim{}_num{}".format(dim, num1d)
     return tuple(constructor_quad_iname(name, d, quadrature_points_per_direction()) for d in range(local_dimension()))
 
 
@@ -135,37 +136,63 @@ def quadrature_weight():
                        dtype=np.float64,
                        shape=(num1d,) * dim,
                        classtag="operator",
-                       dim_tags=",".join(["c"] * dim),
+                       dim_tags=",".join(["f"] * dim),
                        managed=True,
                        potentially_vectorized=True,
                        )
 
     # Precompute it in the constructor
-    instruction(assignee=prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in constructor_quadrature_inames(name))),
-                expression=prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights()), (prim.Variable(i),)) for i in constructor_quadrature_inames(name))),
-                within_inames=frozenset(constructor_quadrature_inames(name)),
+    inames = constructor_quadrature_inames(dim, num1d)
+    instruction(assignee=prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames)),
+                expression=prim.Product(tuple(Subscript(Variable(name_oned_quadrature_weights()), (prim.Variable(i),)) for i in inames)),
+                within_inames=frozenset(inames),
                 kernel="operator",
                 )
 
     return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames()))
 
 
-def define_quadrature_position(name):
-    for i in range(local_dimension()):
-        instruction(expression=Subscript(Variable(name_oned_quadrature_points()), (Variable(quadrature_inames()[i]),)),
-                    assignee=Subscript(Variable(name), (i,)),
-                    forced_iname_deps=frozenset(quadrature_inames()),
-                    forced_iname_deps_is_final=True,
-                    tags=frozenset({"quad"}),
-                    )
+def define_quadrature_position(name, index):
+    instruction(expression=Subscript(Variable(name_oned_quadrature_points()), (Variable(quadrature_inames()[index]),)),
+                assignee=Subscript(Variable(name), (index,)),
+                forced_iname_deps=frozenset(quadrature_inames()),
+                forced_iname_deps_is_final=True,
+                tags=frozenset({"quad"}),
+                )
 
 
 @backend(interface="quad_pos", name="sumfact")
-def pymbolic_quadrature_position():
-    name = 'pos'
-    temporary_variable(name, shape=(local_dimension(),), shape_impl=("fv",))
-    define_quadrature_position(name)
-    return Variable(name)
+def pymbolic_quadrature_position(index):
+    # Return the non-precomputed version
+    if not get_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,))
+
+    assert isinstance(index, int)
+    dim = local_dimension()
+    num1d = quadrature_points_per_direction()
+    name = "quad_points_dim{}_num{}_dir{}".format(dim, num1d, index)
+
+    loopy_class_member(name,
+                       dtype=np.float64,
+                       shape=(num1d,) * dim,
+                       classtag="operator",
+                       dim_tags=",".join(["f"] * dim),
+                       managed=True,
+                       potentially_vectorized=True,
+                       )
+
+    # Precompute it in the constructor
+    inames = constructor_quadrature_inames(dim, num1d)
+    instruction(assignee=prim.Subscript(prim.Variable(name), tuple(prim.Variable(i) for i in inames)),
+                expression=Subscript(Variable(name_oned_quadrature_points()), (prim.Variable(inames[index]))),
+                within_inames=frozenset(inames),
+                kernel="operator",
+                )
+
+    return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames()))
 
 
 @backend(interface="qp_in_cell", name="sumfact")