Compare DBC to numerical solution

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

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)]