Merge branch 'feature/dune-opcounter' into 'master'

Feature/dune opcounter

See merge request !70

dune/perftool/common/opcounter.hh

@@ -15,6 +15,11 @@ namespace oc {
 
 namespace Dune {
 
+  template <typename T>
+  struct IsNumber<oc::OpCounter<T>>
+    : public std::integral_constant<bool, IsNumber<T>::value>{
+  };
+
   template<typename T, int n>
   class FieldVector;
 
@@ -88,6 +93,12 @@ namespace oc {
       return iss;
     }
 
+    friend std::istream& operator>>(std::istream& is, OpCounter& f)
+    {
+      is >> f._v;
+      return is;
+    }
+
     F* data()
     {
       return &_v;

dune/perftool/common/timer.hh

@@ -3,7 +3,7 @@
 
 #include <chrono>
 
-# include <dune/perftool/common/opcounter.hh>
+#include <dune/perftool/common/opcounter.hh>
 
 #define HP_TIMER_OPCOUNTER oc::OpCounter<double>
 

python/dune/perftool/loopy/target.py

+import numpy as np
+
 from dune.perftool.generation import post_include
 
 from dune.perftool.generation.loopy import DuneGlobalArg
@@ -6,25 +8,42 @@ from dune.perftool.loopy.vcl import VCLTypeRegistry
 from dune.perftool.generation import (include_file,
                                       retrieve_cache_functions,
                                       )
+from dune.perftool.options import get_option
 
+from loopy.symbolic import Literal
 from loopy.target import (TargetBase,
                           ASTBuilderBase,
                           DummyHostASTBuilder,
                           )
 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
-
 import pymbolic.primitives as prim
+
+import pytools as pt
+
 import cgen
 
-_registry = {'float32': 'float',
-             'int32': 'int',
-             'float64': 'double',
-             'string': 'std::string',
-             'str': 'std::string'}
+
+def _type_to_op_counter_type(name):
+    return "oc::OpCounter<{}>".format(name)
+
+@pt.memoize
+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 +70,17 @@ class DuneExpressionToCExpressionMapper(ExpressionToCExpressionMapper):
         # the CExpressionToCodeMapper as is => operator/.
         return expr
 
+    def map_constant(self, expr, type_context):
+        ret = ExpressionToCExpressionMapper.map_constant(self, expr, type_context)
+        if get_option('opcounter'):
+            if type_context == "f":
+                _type = _type_to_op_counter_type('float')
+                ret = Literal("{}({})".format(_type, ret.s))
+            if type_context == "d":
+                _type = _type_to_op_counter_type('double')
+                ret = Literal("{}({})".format(_type, ret.s))
+        return ret
+
 
 class DuneCExpressionToCodeMapper(CExpressionToCodeMapper):
     def map_remainder(self, expr, enclosing_prec):
@@ -142,7 +172,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/options.py

@@ -40,6 +40,7 @@ def get_form_compiler_arguments():
     parser.add_argument("--constant-transformation-matrix", action="store_true", help="set option if the jacobian of the transformation is constant on a cell")
     parser.add_argument("--ini-file", type=str, help="An inifile to use. A generated driver will be hard-coded to it, a [formcompiler] section will be used as default values to form compiler arguments (use snake case)")
     parser.add_argument("--timer", action="store_true", help="measure times")
+    parser.add_argument("--opcounter", action="store_true", default=False, help="Count operations. Should only be used with yaspgrid. Timer should be set.")
     parser.add_argument("--project-basedir", type=str, help="The base (build) directory of the dune-perftool project")
     # TODO at some point this help description should be updated
     parser.add_argument("--sumfact", action="store_true", help="Use sumfactorization")

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
@@ -560,6 +565,12 @@ def define_dofestimate(name):
         geo_factor = "6"
     gfs = name_gfs(_driver_data['form'].coefficients()[0].ufl_element())
     ini = name_initree()
+
+    # Assure that gfs in initialized
+    formdata = _driver_data['formdata']
+    x = name_vector(formdata)
+    define_vector(x, formdata)
+
     return ["int generic_dof_estimate =  {} * {}.maxLocalSize();".format(geo_factor, gfs),
             "int {} = {}.get<int>(\"istl.number_of_nnz\", generic_dof_estimate);".format(name, ini)]
 
@@ -1107,17 +1118,15 @@ def dune_solve():
         solve = "{}.apply();".format(snp)
 
     if get_option('timer'):
-        # Necessary includes and defines
-        from dune.perftool.generation import pre_include
+        setup_timer()
         from dune.perftool.generation import post_include
-        pre_include("#define ENABLE_HP_TIMERS", filetag="driver")
-        include_file("dune/perftool/common/timer.hh", filetag="driver")
         post_include("HP_DECLARE_TIMER(total);", filetag="driver")
 
         # Print times after solving
         from dune.perftool.generation import get_global_context_value
         formdatas = get_global_context_value("formdatas")
         print_times = []
+        define_exec()
         for formdata in formdatas:
             lop_name = name_localoperator(formdata)
             timestream = name_timing_stream()
@@ -1125,7 +1134,6 @@ def dune_solve():
         solve = ["HP_TIMER_START(total);",
                  "{}".format(solve),
                  "HP_TIMER_STOP(total);",
-                 "char* exec = argv[0];"
                  "DUMP_TIMER(total, {}, true);".format(timestream),
                  ]
         solve.extend(print_times)
@@ -1310,6 +1318,88 @@ def name_test_fail_variable():
     return name
 
 
+@cached
+def setup_timer():
+    assert(get_option('timer'))
+
+    # Necessary includes and defines
+    from dune.perftool.generation import pre_include
+
+    # TODO check that we are using YASP?
+    if get_option('opcounter'):
+        pre_include("#define ENABLE_COUNTER", filetag="driver")
+    pre_include("#define ENABLE_HP_TIMERS", filetag="driver")
+    include_file("dune/perftool/common/timer.hh", filetag="driver")
+
+
+@preamble
+def define_exec():
+    return "char* exec = argv[0];"
+
+
+@preamble
+def evaluate_residual_timer():
+    formdata = _driver_data['formdata']
+    n_go = name_gridoperator(formdata)
+    v = name_vector(formdata)
+    t_v = type_vector(formdata)
+
+    # Write back times
+    setup_timer()
+    from dune.perftool.generation import post_include
+    post_include("HP_DECLARE_TIMER(residual_evaluation);", filetag="driver")
+    timestream = name_timing_stream()
+    define_exec()
+    print_times = []
+    from dune.perftool.generation import get_global_context_value
+    formdatas = get_global_context_value("formdatas")
+    for formdata in formdatas:
+        lop_name = name_localoperator(formdata)
+        print_times.append("{}.dump_timers({}, argv[0], true);".format(lop_name, timestream))
+
+    evaluation = ["{} r({});".format(t_v, v),
+                  "r=0.0;",
+                  "HP_TIMER_START(residual_evaluation);",
+                  "{}.residual({}, r);".format(n_go, v),
+                  "HP_TIMER_STOP(residual_evaluation);",
+                  "DUMP_TIMER(residual_evaluation, {}, true);".format(timestream)]
+    evaluation.extend(print_times)
+
+    return evaluation
+
+
+@preamble
+def assemble_matrix_timer():
+    formdata = _driver_data['formdata']
+    t_go = type_gridoperator(formdata)
+    n_go = name_gridoperator(formdata)
+    v = name_vector(formdata)
+    t_v = type_vector(formdata)
+
+    # Write back times
+    setup_timer()
+    from dune.perftool.generation import post_include
+    post_include("HP_DECLARE_TIMER(matrix_assembly);", filetag="driver")
+    timestream = name_timing_stream()
+    define_exec()
+    print_times = []
+    from dune.perftool.generation import get_global_context_value
+    formdatas = get_global_context_value("formdatas")
+    for formdata in formdatas:
+        lop_name = name_localoperator(formdata)
+        print_times.append("{}.dump_timers({}, argv[0], true);".format(lop_name, timestream))
+
+    assembly = ["using M = typename {}::Traits::Jacobian;".format(t_go),
+                "M m({});".format(n_go),
+                "HP_TIMER_START(matrix_assembly);",
+                "{}.jacobian({},m);".format(n_go, v),
+                "HP_TIMER_STOP(matrix_assembly);",
+                "DUMP_TIMER(matrix_assembly, {}, true);".format(timestream)]
+    assembly.extend(print_times)
+
+    return assembly
+
+
 @preamble
 def print_residual():
     ini = name_initree()
@@ -1535,9 +1625,15 @@ def generate_driver(formdatas, data):
     # The driver module uses a global dictionary for storing necessary data
     set_driver_data(formdatas, data)
 
-    # The vtkoutput is the generating method that triggers all others.
-    # Alternatively, one could use the `solve` method.
-    if is_stationary():
+    # 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()
+    elif is_stationary():
+        # We could also use solve if we are not interested in visualization
         vtkoutput()
     else:
         solve_instationary()

python/dune/perftool/pdelab/localoperator.py

@@ -9,6 +9,7 @@ from dune.perftool.generation import (backend,
                                       constructor_parameter,
                                       domain,
                                       dump_accumulate_timer,
+                                      generator_factory,
                                       get_backend,
                                       get_global_context_value,
                                       global_context,
@@ -353,7 +354,7 @@ def boundary_predicates(expr, measure, subdomain_id):
         from ufl.checks import is_cellwise_constant
         cellwise_constant = is_cellwise_constant(expr)
         from dune.perftool.pdelab.parameter import intersection_parameter_function
-        intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int')
+        intersection_parameter_function(name, subdomain_data, cellwise_constant, t='int32')
 
         predicates = predicates.union(['{} == {}'.format(name, subdomain_id)])
 
@@ -573,18 +574,27 @@ def name_time_dumper_exec():
     return "exec"
 
 
+@generator_factory(item_tags=("cached",), cache_key_generator=lambda **kw: None)
+def name_example_kernel(name=None):
+    return name
+
+
 class TimerMethod(ClassMember):
     def __init__(self):
         os = name_time_dumper_os()
         reset = name_time_dumper_reset()
         t = name_time_dumper_t()
         ex = name_time_dumper_exec()
+        knl = name_example_kernel()
+        assert(knl is not None)
+
+        # TODO: operator counting only works if alpha_volume_kernel exists
         content = ["template <typename Stream>",
                    "void dump_timers(Stream& {}, char* {}, bool {})".format(os, ex, reset),
                    "{",
                    "  double {} = 0.0;".format(t),
-                   "#ifdef ENABLE_COUNTERS",
-                   "  auto counter = HP_TIMER_OPCOUNTERS({})",
+                   "#ifdef ENABLE_COUNTER",
+                   "  auto counter = HP_TIMER_OPCOUNTERS({});".format(knl),
                    "  counter.reset();",
                    "#endif",
                    ""]
@@ -619,6 +629,7 @@ class LoopyKernelMethod(ClassMember):
             # Start timer
             if add_timings and get_option('timer'):
                 timer_name = assembler_routine_name() + '_kernel'
+                name_example_kernel(name=timer_name)
                 post_include('HP_DECLARE_TIMER({});'.format(timer_name), filetag='operatorfile')
                 content.append('  ' + 'HP_TIMER_START({});'.format(timer_name))
                 dump_accumulate_timer(timer_name)

python/dune/perftool/pdelab/parameter.py

@@ -211,9 +211,11 @@ 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='float64'):
     shape = expr.ufl_element().value_shape()
     shape_impl = ('fv',) * len(shape)
+    from dune.perftool.loopy.target import numpy_to_cpp_dtype
+    t = numpy_to_cpp_dtype(t)
     define_parameter_function_class_member(name, expr, t, shape, True)
     if cellwise_constant:
         evaluate_cellwise_constant_parameter_function(name, restriction)
@@ -223,9 +225,11 @@ 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='float64'):
     shape = expr.ufl_element().value_shape()
     shape_impl = ('fv',) * len(shape)
+    from dune.perftool.loopy.target import numpy_to_cpp_dtype
+    t = numpy_to_cpp_dtype(t)
     define_parameter_function_class_member(name, expr, t, shape, False)
     if cellwise_constant:
         evaluate_intersectionwise_constant_parameter_function(name)

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/opcount_poisson_dg_symdiff.mini

+__name = opcount_poisson_dg_symdiff
+
+cells = 32 32
+extension = 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+extension = vtu
+
+[formcompiler]
+exact_solution_expression = g
+compare_dofs = 1e-1
+compare_l2errorsquared = 1e-4
+opcounter = 1
\ No newline at end of file

test/poisson/poisson_dg_quadrilateral.ufl

+cell = "quadrilateral"
+
+f = Expression("return -2.0*x.size();", cell=cell)
+g = Expression("return x.two_norm2();", on_intersection=True, cell=cell)
+
+V = FiniteElement("DG", cell, 1)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+gamma = 1.0
+theta = 1.0
+
+r = inner(grad(u), grad(v))*dx \
+  + inner(n, avg(grad(u)))*jump(v)*dS \
+  + gamma*jump(u)*jump(v)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - inner(n, grad(u))*v*ds \
+  + gamma*u*v*ds \
+  - theta*u*inner(grad(v), n)*ds \
+  - f*v*dx \
+  + theta*g*inner(grad(v), n)*ds \
+  - gamma*g*v*ds
+
+forms = [r]

test/sumfact/poisson/CMakeLists.txt

@@ -22,6 +22,12 @@ dune_add_formcompiler_system_test(UFLFILE poisson_3d_order2.ufl
                                   )
 
 
+# Operator counting test
+dune_add_formcompiler_system_test(UFLFILE opcount_poisson_2d_order2.ufl
+                                  BASENAME opcount_sumfact_poisson_2d_order2
+                                  INIFILE opcount_poisson_2d_order2.mini
+                                  )
+
 # # 2. Poisson Test Case: DG, f + pure dirichlet
 # dune_add_formcompiler_system_test(UFLFILE poisson_dg.ufl
 #                                  BASENAME sumfact_poisson_dg

test/sumfact/poisson/opcount_poisson_2d_order2.mini

+__name = opcount_sumfact_poisson_2d_order2_{__exec_suffix}
+__exec_suffix = {diff_suffix}
+
+diff_suffix = symdiff
+
+cells = 8 8
+extension = 1. 1.
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_ref
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 0
+exact_solution_expression = g
+compare_l2errorsquared = 1e-8
+sumfact = 1
+opcounter = 1

test/sumfact/poisson/opcount_poisson_2d_order2.ufl

+cell = "quadrilateral"
+
+f = Expression("Dune::FieldVector<oc::OpCounter<double>,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*exp(-1.*c.two_norm2());", cell=cell)
+g = Expression("Dune::FieldVector<oc::OpCounter<double>,2> c(0.5); c-= x; return exp(-1.*c.two_norm2());", cell=cell)
+
+V = FiniteElement("CG", cell, 2, dirichlet_expression=g)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+forms = [(inner(grad(u), grad(v)) - f*v)*dx]