From fd073fedf407a901a63a4b33969e5b8b2b304af9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Thu, 11 Aug 2016 09:40:44 +0200
Subject: [PATCH] Make it possible to compare l2errorsquared

---
 python/dune/perftool/options.py       | 14 ++++++-
 python/dune/perftool/pdelab/driver.py | 57 +++++++++++++++++++++++----
 test/nonlinear/CMakeLists.txt         |  2 +-
 test/nonlinear/nonlinear_dg.mini      | 11 ++++++
 4 files changed, 73 insertions(+), 11 deletions(-)
 create mode 100644 test/nonlinear/nonlinear_dg.mini

diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index e4b51c3a..ff563cfb 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -15,8 +15,18 @@ def get_form_compiler_arguments():
     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(
         "--exact-solution-expression",
-        nargs = 2,
-        help="requires two arguments: the name of the expression corresponding to solution in the ufl file and the maximal allowed error (between difference of the dof vectors of num. solution and interpolated exact solution")
+        type=str,
+        help="name of the exact solution expression in the ufl file")
+    parser.add_argument(
+        "--compare-dofs",
+        type=str,
+        help="maximal allowed maximum error of difference between degrees of freedom vectors of numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)"
+    )
+    parser.add_argument(
+        "--compare-l2errorsquared",
+        type=str,
+        help="maximal allowed l2 error squared of difference between numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)"
+    )
     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 12c6b674..f246e0d9 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -749,7 +749,7 @@ def name_boundary_function(expr):
 def name_solution_function(expr):
     # Get expression representing solution
     from dune.perftool.generation import get_global_context_value
-    solution_expression_name = get_option("exact_solution_expression")[0]
+    solution_expression_name = get_option("exact_solution_expression")
     solution_expression = get_global_context_value("data").object_by_name[solution_expression_name]
 
     # Create function just like it is done for boundary_function
@@ -924,19 +924,55 @@ def name_vtkfile():
 
 
 @preamble
-def exact_solution():
+def compare_dofs():
     v = name_vector()
     solution = name_vector_from_solution_expression()
-
-    return ["using Dune::PDELab::Backend::native;",
+    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("exact_solution_expression")[1]),
+            "if (maxerror>{})".format(get_option("compare_dofs")),
+            "  throw;"]
+
+
+@preamble
+def define_discrete_grid_function(gfs, vector_name, dgf_name):
+    return ["typedef Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})> DGF;".format(gfs, vector_name),
+            "DGF {}({},{});".format(dgf_name, gfs, vector_name)]
+
+
+@symbol
+def name_discrete_grid_function(gfs, vector_name):
+    dgf_name = vector_name + "_dgf"
+    define_discrete_grid_function(gfs, vector_name, dgf_name)
+    return dgf_name
+
+
+@preamble
+def compare_L2_squared():
+    element = _form.coefficients()[0].ufl_element()
+    v = name_vector()
+    gfs = name_gfs(element)
+    vdgf = name_discrete_grid_function(gfs, v)
+    solution_function = name_solution_function(element)
+
+    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",
+            "typedef DifferenceSquaredAdapter<decltype({}),decltype({})> DifferenceSquared;".format(solution_function, vdgf),
+            "DifferenceSquared differencesquared({},{});".format(solution_function, vdgf),
+            "typename DifferenceSquared::Traits::RangeType l2errorsquared(0.0);",
+            "Dune::PDELab::integrateGridFunction(differencesquared,l2errorsquared,10);",
+            "std::cout << \"l2errorsquared: \" << l2errorsquared << std::endl;",
+            "if (l2errorsquared>{})".format(get_option("compare_l2errorsquared")),
             "  throw;"]
 
 
@@ -1024,7 +1060,12 @@ def vtkoutput():
         from ufl import MixedElement, VectorElement, TensorElement
         if isinstance(_form.coefficients()[0].ufl_element(),  (MixedElement, VectorElement, TensorElement)):
             raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
-        exact_solution()
+
+        if get_option("compare_dofs"):
+            compare_dofs()
+        if get_option("compare_l2errorsquared"):
+            compare_L2_squared()
+
 
     return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vec, predicate),
             "{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]
diff --git a/test/nonlinear/CMakeLists.txt b/test/nonlinear/CMakeLists.txt
index 5ebedaa8..a78cccd0 100644
--- a/test/nonlinear/CMakeLists.txt
+++ b/test/nonlinear/CMakeLists.txt
@@ -11,7 +11,7 @@ dune_add_system_test(TARGET nonlinear_symdiff
 
 add_generated_executable(UFLFILE nonlinear_dg.ufl
                          TARGET nonlinear_dg
-                         FORM_COMPILER_ARGS --exact-solution-expression g 1e-2
+                         FORM_COMPILER_ARGS --exact-solution-expression g --compare-dofs 1e-1 --compare-l2errorsquared 1e-3
                          )
 
 dune_add_system_test(TARGET nonlinear_dg
diff --git a/test/nonlinear/nonlinear_dg.mini b/test/nonlinear/nonlinear_dg.mini
new file mode 100644
index 00000000..d2b26e27
--- /dev/null
+++ b/test/nonlinear/nonlinear_dg.mini
@@ -0,0 +1,11 @@
+__name = nonlinear_dg
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = nonlinear_dg
+extension = vtu
-- 
GitLab