diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py index 82c53c1da31e890ad39fc321bea6e495b73ecab5..1aeb7f272d52b7c23ac361b9df132fa1381d36d6 100644 --- a/python/dune/perftool/options.py +++ b/python/dune/perftool/options.py @@ -13,6 +13,10 @@ def get_form_compiler_arguments(): parser.add_argument("--operator-file", type=str, help="The filename for the generated local operator header") parser.add_argument('--version', action='version', version='%(prog)s 0.1') parser.add_argument("--numerical-jacobian", action="store_true", help="use numerical jacobians (only makes sense, if uflpdelab for some reason fails to generate analytic jacobians)") + parser.add_argument( + "--dirichlet-is-exact-solution", + action="store_true", + help="test if solution is correct by providing the exact solution as dirichlet boundary condition") parser.add_argument("--interactive", action="store_true", help="whether the optimization process should be guided interactively (also useful for debugging)") # These are the options that I deemed necessary in uflpdelab diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py index 4fe77cdfb46c61e13c6104caa8ef2099fe279d5b..1d01247dc175e137b08ebdc0e279f8497ae3827a 100644 --- a/python/dune/perftool/pdelab/driver.py +++ b/python/dune/perftool/pdelab/driver.py @@ -11,6 +11,7 @@ from dune.perftool.generation import (generator_factory, symbol, preamble, ) +from dune.perftool.options import get_option def has_constraints(element): @@ -756,6 +757,10 @@ def name_vector(): maybe_interpolate_vector("x") return "x" +@symbol +def name_solution_from_dirichlet_vector(): + interpolate_vector("solution") + return "solution" @preamble def typedef_linearsolver(name): @@ -849,6 +854,24 @@ def name_vtkfile(): return "vtkfile" +@preamble +def exact_solution_from_dirichlet(): + v = name_vector() + solution = name_solution_from_dirichlet_vector(); + + return ["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 << maxerror << std::endl;", + "", + "if (maxerror>1e-14)", + " throw;"] + + + @preamble def print_residual(): ini = name_initree() @@ -928,5 +951,8 @@ def vtkoutput(): dune_solve() print_residual() print_matrix() + if get_option("dirichlet_is_exact_solution"): + exact_solution_from_dirichlet() + return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vec, predicate), "{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]