Implement OpCounter in local operator

dune/perftool/common/opcounter.hh

@@ -88,6 +88,12 @@ namespace oc {
       return iss;
     }
 
+    friend std::istream& operator>>(std::istream& is, OpCounter& f)
+    {
+      is >> f._v;
+      return is;
+    }
+
     F* data()
     {
       return &_v;

python/dune/perftool/loopy/target.py

+import numpy as np
+
 from dune.perftool.generation import post_include
 
 from dune.perftool.generation.loopy import DuneGlobalArg
 from dune.perftool.loopy.temporary import DuneTemporaryVariable
 from dune.perftool.loopy.vcl import VCLTypeRegistry
+from dune.perftool.options import get_option
 from dune.perftool.generation import (include_file,
                                       retrieve_cache_functions,
                                       )
@@ -13,6 +16,7 @@ from loopy.target import (TargetBase,
                           )
 from loopy.target.c import CASTBuilder
 from loopy.target.c.codegen.expression import ExpressionToCExpressionMapper, CExpressionToCodeMapper
+from loopy.tools import is_integer
 from loopy.types import NumpyType
 
 from pymbolic.mapper.stringifier import PREC_NONE
@@ -20,11 +24,23 @@ from pymbolic.mapper.stringifier import PREC_NONE
 import pymbolic.primitives as prim
 import cgen
 
-_registry = {'float32': 'float',
-             'int32': 'int',
-             'float64': 'double',
-             'string': 'std::string',
-             'str': 'std::string'}
+
+def _type_to_op_counter_type(type):
+    return "oc::OpCounter<{}>".format(type)
+
+
+def numpy_to_cpp_dtype(key):
+    _registry = {'float32': 'float',
+                 'int32': 'int',
+                 'float64': 'double',
+                 'string': 'std::string',
+                 'str': 'std::string'}
+
+    if get_option('opcounter'):
+        _registry['float32'] = _type_to_op_counter_type('float')
+        _registry['float64'] = _type_to_op_counter_type('double')
+
+    return _registry[key]
 
 
 class DuneExpressionToCExpressionMapper(ExpressionToCExpressionMapper):
@@ -51,6 +67,31 @@ class DuneExpressionToCExpressionMapper(ExpressionToCExpressionMapper):
         # the CExpressionToCodeMapper as is => operator/.
         return expr
 
+    def map_constant(self, expr, type_context):
+        if isinstance(expr, (complex, np.complexfloating)):
+            return ExpressionToCExpressionMapper.map_constant(self, expr, type_context)
+
+        else:
+            from loopy.symbolic import Literal
+            if type_context == "f":
+                if get_option('opcounter'):
+                    type = _type_to_op_counter_type('float')
+                    return Literal("{}({})".format(type,repr(float(expr))+'f'))
+                return Literal(repr(float(expr))+"f")
+            elif type_context == "d":
+                if get_option('opcounter'):
+                    type = _type_to_op_counter_type('double')
+                    return Literal("{}({})".format(type,repr(float(expr))))
+                return Literal(repr(float(expr)))
+            elif type_context == "i":
+                return int(expr)
+            else:
+                if is_integer(expr):
+                    return int(expr)
+
+                raise RuntimeError("don't know how to generate code "
+                        "for constant '%s'" % expr)
+
 
 class DuneCExpressionToCodeMapper(CExpressionToCodeMapper):
     def map_remainder(self, expr, enclosing_prec):
@@ -142,7 +183,7 @@ class DuneTarget(TargetBase):
             include_file("dune/perftool/vectorclass/vectorclass.h", filetag="operatorfile")
             return VCLTypeRegistry.names[dtype.dtype]
         else:
-            return _registry[dtype.dtype.name]
+            return numpy_to_cpp_dtype(dtype.dtype.name)
 
     def is_vector_dtype(self, dtype):
         return False

python/dune/perftool/loopy/temporary.py

@@ -9,12 +9,13 @@ import numpy
 
 
 def _temporary_type(shape_impl, shape, first=True):
+    from dune.perftool.loopy.target import numpy_to_cpp_dtype
     if len(shape_impl) == 0:
-        return 'double'
+        return numpy_to_cpp_dtype('float64')
     if shape_impl[0] == 'arr':
         if not first or len(set(shape_impl)) != 1:
             raise PerftoolLoopyError("We do not allow mixing of C++ containers and plain C arrays, for reasons of mental sanity")
-        return 'double'
+        return numpy_to_cpp_dtype('float64')
     if shape_impl[0] == 'vec':
         return "std::vector<{}>".format(_temporary_type(shape_impl[1:], shape[1:], first=False))
     if shape_impl[0] == 'fv':
@@ -22,7 +23,8 @@ def _temporary_type(shape_impl, shape, first=True):
     if shape_impl[0] == 'fm':
         # For now, no field matrices of weird stuff...
         assert len(shape) == 2
-        return "Dune::FieldMatrix<double, {}, {}>".format(shape[0], shape[1])
+        type = numpy_to_cpp_dtype('float64')
+        return "Dune::FieldMatrix<{}, {}, {}>".format(type, shape[0], shape[1])
 
 
 def default_declaration(name, shape=(), shape_impl=()):

python/dune/perftool/pdelab/driver.py

@@ -181,7 +181,8 @@ def typedef_grid(name):
         set_option('diagonal_transformation_matrix', True)
         set_option('constant_transformation_matrix', True)
 
-        gridt = "Dune::YaspGrid<{}>".format(dim)
+        range_type = type_range()
+        gridt = "Dune::YaspGrid<{}, Dune::EquidistantCoordinates<{},dim>>".format(dim, range_type)
         include_file("dune/grid/yaspgrid.hh", filetag="driver")
     else:
         if any(_driver_data['form'].ufl_cell().cellname() in x for x in ["triangle", "tetrahedron"]):
@@ -285,12 +286,16 @@ def type_domainfield():
 
 @preamble
 def typedef_range(name):
-    return "using {} = double;".format(name)
+    if get_option('opcounter'):
+        setup_timer()
+        return "using {} = oc::OpCounter<double>;".format(name)
+    else:
+        return "using {} = double;".format(name)
 
 
 def type_range():
-    typedef_range("R")
-    return "R"
+    typedef_range("RangeType")
+    return "RangeType"
 
 
 @preamble
@@ -1352,8 +1357,7 @@ def evaluate_residual_timer():
         lop_name = name_localoperator(formdata)
         print_times.append("{}.dump_timers({}, argv[0], true);".format(lop_name, timestream))
 
-    evaluation = ["using Dune::PDELab::Backend::native;",
-                  "{} r({});".format(t_v, v),
+    evaluation = ["{} r({});".format(t_v, v),
                   "r=0.0;",
                   "HP_TIMER_START(residual_evaluation);",
                   "{}.residual({}, r);".format(n_go, v),
@@ -1623,6 +1627,8 @@ def generate_driver(formdatas, data):
 
     # Entrypoint for driver generation
     if get_option("opcounter"):
+        assert(any(_driver_data['form'].ufl_cell().cellname() in x for x in
+                   ["vertex", "interval", "quadrilateral", "hexahedron"]))
         # In case of operator conunting we only assemble the matrix and evaluate the residual
         assemble_matrix_timer()
         evaluate_residual_timer()

python/dune/perftool/pdelab/parameter.py

@@ -10,6 +10,7 @@ from dune.perftool.generation import (class_basename,
                                       preamble,
                                       temporary_variable
                                       )
+from dune.perftool.loopy.target import numpy_to_cpp_dtype
 from dune.perftool.pdelab.geometry import (name_cell,
                                            name_dimension,
                                            name_intersection,
@@ -211,7 +212,7 @@ def construct_nested_fieldvector(t, shape):
 
 
 @kernel_cached
-def cell_parameter_function(name, expr, restriction, cellwise_constant, t='double'):
+def cell_parameter_function(name, expr, restriction, cellwise_constant, t=numpy_to_cpp_dtype('float64')):
     shape = expr.ufl_element().value_shape()
     shape_impl = ('fv',) * len(shape)
     define_parameter_function_class_member(name, expr, t, shape, True)
@@ -223,7 +224,7 @@ def cell_parameter_function(name, expr, restriction, cellwise_constant, t='doubl
 
 
 @kernel_cached
-def intersection_parameter_function(name, expr, cellwise_constant, t='double'):
+def intersection_parameter_function(name, expr, cellwise_constant, t=numpy_to_cpp_dtype('float64')):
     shape = expr.ufl_element().value_shape()
     shape_impl = ('fv',) * len(shape)
     define_parameter_function_class_member(name, expr, t, shape, False)

test/poisson/CMakeLists.txt

@@ -48,7 +48,14 @@ dune_add_formcompiler_system_test(UFLFILE poisson_cellwise_constant.ufl
                                   INIFILE poisson_cellwise_constant.mini
                                   )
 
-# Add the reference code with the PDELab localoperator that produced
+# 8. Poisson with operator counting
+dune_add_formcompiler_system_test(UFLFILE poisson_dg_quadrilateral.ufl
+                                  BASENAME opcount_poisson_dg_symdiff
+                                  INIFILE opcount_poisson_dg_symdiff.mini
+                                  )
+
+
+
 # the reference vtk file
 add_executable(poisson_dg_ref reference_main.cc)
 set_target_properties(poisson_dg_ref PROPERTIES EXCLUDE_FROM_ALL 1)

test/poisson/dimension-grid-variations/poisson_2d_dg_quadrilateral.mini

@@ -11,4 +11,3 @@ extension = vtu
 exact_solution_expression = g
 compare_dofs = 1e-1
 compare_l2errorsquared = 1e-4
-opcounter = 1
\ No newline at end of file