diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py index abce3936ded3e07e73971c0ca67e5c4fefb7c5e9..b12a2fd90e2c4cb7491d0ad2f84a62069bcdf183 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 0000000000000000000000000000000000000000..c02147a94655a5c9bba0716d4e98b2cd7d84e019 --- /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)]