From 5d43bb6729c58ae62ebf4a1ba115f8398c93c1c6 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Wed, 21 Jun 2017 13:28:49 +0200
Subject: [PATCH] Extract any error calculation from driver module
---
.../dune/perftool/pdelab/driver/__init__.py | 214 +--------------
python/dune/perftool/pdelab/driver/error.py | 244 ++++++++++++++++++
2 files changed, 246 insertions(+), 212 deletions(-)
create mode 100644 python/dune/perftool/pdelab/driver/error.py
diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py
index abce3936..b12a2fd9 100644
--- a/python/dune/perftool/pdelab/driver/__init__.py
+++ b/python/dune/perftool/pdelab/driver/__init__.py
@@ -890,35 +890,6 @@ def name_boundary_function(expr):
return _name_boundary_function(element, boundary)
-@cached
-def name_solution_function(tree_path=()):
- from dune.perftool.generation import get_global_context_value
- expr = get_global_context_value("data").object_by_name[get_option("exact_solution_expression")]
-
- from dune.perftool.ufl.execution import Expression, split_expression
- from ufl.classes import ListTensor
- for i in tree_path:
- if isinstance(expr, Expression):
- expr = split_expression(expr)[int(i)]
- elif isinstance(expr, tuple):
- expr = expr[int(i)]
- elif isinstance(expr, ListTensor):
- expr = expr.ufl_operands[int(i)]
- else:
- raise TypeError("No idea how to split {}".format(type(expr)))
-
- from dune.perftool.generation import get_global_context_value
- try:
- name = get_global_context_value("data").object_names[id(expr)]
- except KeyError:
- from dune.perftool.generation import get_counter
- name = 'generic_{}'.format(get_counter('__generic'))
-
- define_boundary_function(expr, name)
-
- return name
-
-
@preamble
def interpolate_vector(name, formdata):
define_vector(name, formdata)
@@ -931,19 +902,6 @@ def interpolate_vector(name, formdata):
)
-@preamble
-def interpolate_solution_expression(name):
- formdata = _driver_data['formdata']
- define_vector(name, formdata)
- element = get_trial_element()
- sol = name_solution_function()
- gfs = name_gfs(element)
- return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
- gfs,
- name,
- )
-
-
def maybe_interpolate_vector(name, formdata):
element = get_trial_element()
if has_constraints(element):
@@ -958,11 +916,6 @@ def name_vector(formdata):
return name
-def name_vector_from_solution_expression():
- interpolate_solution_expression("solution")
- return "solution"
-
-
@preamble
def typedef_linearsolver(name):
include_file("dune/pdelab/backend/istl.hh", filetag="driver")
@@ -1250,171 +1203,6 @@ def name_vtkfile():
return "vtkfile"
-@preamble
-def compare_dofs():
- v = name_vector(_driver_data['formdata'])
- solution = name_vector_from_solution_expression()
- fail = name_test_fail_variable()
- return ["",
- "// Maximum norm of difference between dof vectors of the",
- "// numerical solution and the interpolation of the exact solution",
- "using Dune::PDELab::Backend::native;",
- "double maxerror(0.0);",
- "for (std::size_t i=0; i<native({}).size(); ++i)".format(v),
- " if (std::abs(native({})[i]-native({})[i]) > maxerror)".format(v, solution),
- " maxerror = std::abs(native({})[i]-native({})[i]);".format(v, solution),
- "std::cout << std::endl << \"Maxerror: \" << maxerror << std::endl;",
- "if (maxerror>{})".format(get_option("compare_dofs")),
- " {} = true;".format(fail)]
-
-
-def type_discrete_grid_function(gfs):
- return "DGF_{}".format(gfs.upper())
-
-
-@preamble
-def define_discrete_grid_function(gfs, vector_name, dgf_name):
- dgf_type = type_discrete_grid_function(gfs)
- return ["using {} = Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})>;".format(dgf_type, gfs, vector_name),
- "{} {}({},{});".format(dgf_type, dgf_name, gfs, vector_name)]
-
-
-def name_discrete_grid_function(gfs, vector_name):
- dgf_name = gfs + "_dgf"
- define_discrete_grid_function(gfs, vector_name, dgf_name)
- return dgf_name
-
-
-def type_subgfs(element, tree_path):
- include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
- gfs = type_gfs(element)
- return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(tree_path))
-
-
-@preamble
-def define_subgfs(name, element, tree_path):
- t = type_subgfs(element, tree_path)
- gfs = name_gfs(element)
-
- return '{} {}({});'.format(t, name, gfs)
-
-
-def name_subgfs(element, tree_path):
- gfs = name_gfs(element)
- if len(tree_path) == 0:
- return gfs
- else:
- name = "subgfs_{}_{}".format(gfs, '_'.join(tree_path))
- define_subgfs(name, element, tree_path)
- return name
-
-
-@preamble
-def typedef_difference_squared_adapter(name, tree_path):
- element = get_trial_element()
- formdata = _driver_data['formdata']
-
- solution_function = name_solution_function(tree_path)
- gfs = name_subgfs(element, tree_path)
- vector = name_vector(formdata)
- dgf = name_discrete_grid_function(gfs, vector)
-
- return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, solution_function, dgf)
-
-
-def type_difference_squared_adapter(tree_path):
- t = 'DifferenceSquared_{}'.format('_'.join(tree_path))
- typedef_difference_squared_adapter(t, tree_path)
- return t
-
-
-@preamble
-def define_difference_squared_adapter(name, tree_path):
- element = get_trial_element()
- formdata = _driver_data['formdata']
-
- t = type_difference_squared_adapter(tree_path)
- solution_function = name_solution_function(tree_path)
- gfs = name_subgfs(element, tree_path)
- vector = name_vector(formdata)
- dgf = name_discrete_grid_function(gfs, vector)
-
- return '{} {}({}, {});'.format(t, name, solution_function, dgf)
-
-
-def name_difference_squared_adapter(tree_path):
- name = 'differencesquared_{}'.format('_'.join(tree_path))
- define_difference_squared_adapter(name, tree_path)
- return name
-
-
-@preamble
-def _accumulate_L2_squared(tree_path):
- dsa = name_difference_squared_adapter(tree_path)
- t_dsa = type_difference_squared_adapter(tree_path)
- accum_error = name_accumulated_L2_error()
-
- include_file("dune/pdelab/gridfunctionspace/gridfunctionadapter.hh", filetag="driver")
- include_file("dune/pdelab/common/functionutilities.hh", filetag="driver")
-
- return ["{",
- " // L2 error squared of difference between numerical",
- " // solution and the interpolation of exact solution",
- " // only subgfs with treepath: {}".format(', '.join(tree_path)),
- " typename {}::Traits::RangeType err(0.0);".format(t_dsa),
- " Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
- " {} += err;".format(accum_error),
- " std::cout << \"Accumlated L2 Error for tree path {}: \" << err << std::endl;".format(tree_path),
- "}"]
-
-
-def accumulate_L2_squared(tree_path=()):
- element = get_trial_element()
- from ufl.functionview import select_subelement
- from ufl.classes import MultiIndex, FixedIndex
- element = select_subelement(element, MultiIndex(tuple(FixedIndex(int(i)) for i in tree_path)))
-
- from ufl import MixedElement
- if isinstance(element, MixedElement):
- for i, subel in enumerate(element.sub_elements()):
- accumulate_L2_squared(tree_path + (str(i),))
- else:
- _accumulate_L2_squared(tree_path)
-
-
-@preamble
-def define_accumulated_L2_error(name):
- return "double {}(0.0);".format(name)
-
-
-def name_accumulated_L2_error():
- name = 'l2error'
- define_accumulated_L2_error(name)
- return name
-
-
-@preamble
-def compare_L2_squared():
- accumulate_L2_squared()
-
- accum_error = name_accumulated_L2_error()
- fail = name_test_fail_variable()
- return ["std::cout << \"l2errorsquared: \" << {} << std::endl;".format(accum_error),
- "if (std::isnan({0}) or std::abs({0})>{1})".format(accum_error, get_option("compare_l2errorsquared")),
- " {} = true;".format(fail)]
-
-
-@preamble
-def define_test_fail_variable(name):
- return 'bool {}(false);'.format(name)
-
-
-def name_test_fail_variable():
- name = "testfail"
- define_test_fail_variable(name)
- return name
-
-
@preamble
def print_residual():
ini = name_initree()
@@ -1502,6 +1290,7 @@ def vtkoutput():
print_residual()
print_matrix()
+ from dune.perftool.pdelab.driver.error import compare_dofs, compare_L2_squared
if get_option("exact_solution_expression"):
if get_option("compare_dofs"):
compare_dofs()
@@ -1639,6 +1428,7 @@ def solve_instationary():
@preamble
def return_statement():
+ from dune.perftool.pdelab.driver.error import name_test_fail_variable
fail = name_test_fail_variable()
return "return {};".format(fail)
diff --git a/python/dune/perftool/pdelab/driver/error.py b/python/dune/perftool/pdelab/driver/error.py
new file mode 100644
index 00000000..c02147a9
--- /dev/null
+++ b/python/dune/perftool/pdelab/driver/error.py
@@ -0,0 +1,244 @@
+""" Generator functions to calculate errors in the driver """
+
+from dune.perftool.generation import (cached,
+ include_file,
+ preamble,
+ )
+from dune.perftool.options import get_option
+
+
+@preamble
+def define_test_fail_variable(name):
+ return 'bool {}(false);'.format(name)
+
+
+def name_test_fail_variable():
+ name = "testfail"
+ define_test_fail_variable(name)
+ return name
+
+
+def name_vector_from_solution_expression():
+ interpolate_solution_expression("solution")
+ return "solution"
+
+
+@preamble
+def interpolate_solution_expression(name):
+ from dune.perftool.pdelab.driver import (define_vector,
+ get_formdata,
+ get_trial_element,
+ name_gfs,
+ )
+ formdata = get_formdata()
+ define_vector(name, formdata)
+ element = get_trial_element()
+ sol = name_solution_function()
+ gfs = name_gfs(element)
+ return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
+ gfs,
+ name,
+ )
+
+
+@preamble
+def compare_dofs():
+ from dune.perftool.pdelab.driver import (get_formdata,
+ name_vector,
+ )
+
+ v = name_vector(get_formdata())
+ solution = name_vector_from_solution_expression()
+ fail = name_test_fail_variable()
+
+ return ["",
+ "// Maximum norm of difference between dof vectors of the",
+ "// numerical solution and the interpolation of the exact solution",
+ "using Dune::PDELab::Backend::native;",
+ "double maxerror(0.0);",
+ "for (std::size_t i=0; i<native({}).size(); ++i)".format(v),
+ " if (std::abs(native({})[i]-native({})[i]) > maxerror)".format(v, solution),
+ " maxerror = std::abs(native({})[i]-native({})[i]);".format(v, solution),
+ "std::cout << std::endl << \"Maxerror: \" << maxerror << std::endl;",
+ "if (maxerror>{})".format(get_option("compare_dofs")),
+ " {} = true;".format(fail)]
+
+
+
+def type_discrete_grid_function(gfs):
+ return "DGF_{}".format(gfs.upper())
+
+
+@preamble
+def define_discrete_grid_function(gfs, vector_name, dgf_name):
+ dgf_type = type_discrete_grid_function(gfs)
+ return ["using {} = Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})>;".format(dgf_type, gfs, vector_name),
+ "{} {}({},{});".format(dgf_type, dgf_name, gfs, vector_name)]
+
+
+def name_discrete_grid_function(gfs, vector_name):
+ dgf_name = gfs + "_dgf"
+ define_discrete_grid_function(gfs, vector_name, dgf_name)
+ return dgf_name
+
+
+def type_subgfs(element, tree_path):
+ from dune.perftool.pdelab.driver import type_gfs
+ include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
+ gfs = type_gfs(element)
+ return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(tree_path))
+
+
+@preamble
+def define_subgfs(name, element, tree_path):
+ t = type_subgfs(element, tree_path)
+
+ from dune.perftool.pdelab.driver import name_gfs
+ gfs = name_gfs(element)
+
+ return '{} {}({});'.format(t, name, gfs)
+
+
+def name_subgfs(element, tree_path):
+ from dune.perftool.pdelab.driver import name_gfs
+ gfs = name_gfs(element)
+ if len(tree_path) == 0:
+ return gfs
+ else:
+ name = "subgfs_{}_{}".format(gfs, '_'.join(tree_path))
+ define_subgfs(name, element, tree_path)
+ return name
+
+
+@cached
+def name_solution_function(tree_path=()):
+ from dune.perftool.generation import get_global_context_value
+ expr = get_global_context_value("data").object_by_name[get_option("exact_solution_expression")]
+
+ from dune.perftool.ufl.execution import Expression, split_expression
+ from ufl.classes import ListTensor
+ for i in tree_path:
+ if isinstance(expr, Expression):
+ expr = split_expression(expr)[int(i)]
+ elif isinstance(expr, tuple):
+ expr = expr[int(i)]
+ elif isinstance(expr, ListTensor):
+ expr = expr.ufl_operands[int(i)]
+ else:
+ raise TypeError("No idea how to split {}".format(type(expr)))
+
+ from dune.perftool.generation import get_global_context_value
+ try:
+ name = get_global_context_value("data").object_names[id(expr)]
+ except KeyError:
+ from dune.perftool.generation import get_counter
+ name = 'generic_{}'.format(get_counter('__generic'))
+
+ from dune.perftool.pdelab.driver import define_boundary_function
+ define_boundary_function(expr, name)
+
+ return name
+
+
+@preamble
+def typedef_difference_squared_adapter(name, tree_path):
+ from dune.perftool.pdelab.driver import (get_formdata,
+ get_trial_element,
+ name_vector,
+ )
+ element = get_trial_element()
+ formdata = get_formdata()
+
+ solution_function = name_solution_function(tree_path)
+ gfs = name_subgfs(element, tree_path)
+ vector = name_vector(formdata)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return 'using {} = Dune::PDELab::DifferenceSquaredAdapter<decltype({}), decltype({})>;'.format(name, solution_function, dgf)
+
+
+def type_difference_squared_adapter(tree_path):
+ t = 'DifferenceSquared_{}'.format('_'.join(tree_path))
+ typedef_difference_squared_adapter(t, tree_path)
+ return t
+
+
+@preamble
+def define_difference_squared_adapter(name, tree_path):
+ from dune.perftool.pdelab.driver import (get_formdata,
+ get_trial_element,
+ name_vector,
+ )
+ element = get_trial_element()
+ formdata = get_formdata()
+
+ t = type_difference_squared_adapter(tree_path)
+ solution_function = name_solution_function(tree_path)
+ gfs = name_subgfs(element, tree_path)
+ vector = name_vector(formdata)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return '{} {}({}, {});'.format(t, name, solution_function, dgf)
+
+
+def name_difference_squared_adapter(tree_path):
+ name = 'differencesquared_{}'.format('_'.join(tree_path))
+ define_difference_squared_adapter(name, tree_path)
+ return name
+
+
+@preamble
+def _accumulate_L2_squared(tree_path):
+ dsa = name_difference_squared_adapter(tree_path)
+ t_dsa = type_difference_squared_adapter(tree_path)
+ accum_error = name_accumulated_L2_error()
+
+ include_file("dune/pdelab/gridfunctionspace/gridfunctionadapter.hh", filetag="driver")
+ include_file("dune/pdelab/common/functionutilities.hh", filetag="driver")
+
+ return ["{",
+ " // L2 error squared of difference between numerical",
+ " // solution and the interpolation of exact solution",
+ " // only subgfs with treepath: {}".format(', '.join(tree_path)),
+ " typename {}::Traits::RangeType err(0.0);".format(t_dsa),
+ " Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
+ " {} += err;".format(accum_error),
+ " std::cout << \"Accumlated L2 Error for tree path {}: \" << err << std::endl;".format(tree_path),
+ "}"]
+
+
+def accumulate_L2_squared(tree_path=()):
+ from dune.perftool.pdelab.driver import get_trial_element
+ element = get_trial_element()
+ from ufl.functionview import select_subelement
+ from ufl.classes import MultiIndex, FixedIndex
+ element = select_subelement(element, MultiIndex(tuple(FixedIndex(int(i)) for i in tree_path)))
+
+ from ufl import MixedElement
+ if isinstance(element, MixedElement):
+ for i, subel in enumerate(element.sub_elements()):
+ accumulate_L2_squared(tree_path + (str(i),))
+ else:
+ _accumulate_L2_squared(tree_path)
+
+
+@preamble
+def define_accumulated_L2_error(name):
+ return "double {}(0.0);".format(name)
+
+
+def name_accumulated_L2_error():
+ name = 'l2error'
+ define_accumulated_L2_error(name)
+ return name
+
+
+@preamble
+def compare_L2_squared():
+ accumulate_L2_squared()
+
+ accum_error = name_accumulated_L2_error()
+ fail = name_test_fail_variable()
+ return ["std::cout << \"l2errorsquared: \" << {} << std::endl;".format(accum_error),
+ "if (std::isnan({0}) or std::abs({0})>{1})".format(accum_error, get_option("compare_l2errorsquared")),
+ " {} = true;".format(fail)]
--
GitLab