Formcompiler argument marking exact solution

NOTE: Works only for scalar equations

python/dune/perftool/options.py

@@ -14,9 +14,9 @@ def get_form_compiler_arguments():
     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")
+        "--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")
     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

@@ -706,6 +706,7 @@ def name_boundary_lambda(boundary, name):
 def define_boundary_function(boundary, name):
     gv = name_leafview()
     lambdaname = name_boundary_lambda(boundary, name)
+    include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
     return "auto {} = Dune::PDELab::makeGridFunctionFromCallable({}, {});".format(name,
                                                                                   gv,
                                                                                   lambdaname,
@@ -744,6 +745,20 @@ def name_boundary_function(expr):
     return name
 
 
+@symbol
+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 = get_global_context_value("data").object_by_name[solution_expression_name]
+
+    # Create function just like it is done for boundary_function
+    name = "solution_function"
+    define_boundary_function(solution_expression, name)
+
+    return name
+
+
 @preamble
 def interpolate_vector(name):
     define_vector(name)
@@ -756,6 +771,18 @@ def interpolate_vector(name):
                                                            )
 
 
+@preamble
+def interpolate_solution_expression(name):
+    define_vector(name)
+    element = _form.coefficients()[0].ufl_element()
+    sol = name_solution_function(element)
+    gfs = name_gfs(element)
+    return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
+                                                           gfs,
+                                                           name,
+                                                           )
+
+
 def maybe_interpolate_vector(name):
     element = _form.coefficients()[0].ufl_element()
     if has_constraints(element):
@@ -769,11 +796,13 @@ def name_vector():
     maybe_interpolate_vector("x")
     return "x"
 
+
 @symbol
-def name_solution_from_dirichlet_vector():
-    interpolate_vector("solution")
+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")
@@ -895,9 +924,9 @@ def name_vtkfile():
 
 
 @preamble
-def exact_solution_from_dirichlet():
+def exact_solution():
     v = name_vector()
-    solution = name_solution_from_dirichlet_vector();
+    solution = name_vector_from_solution_expression()
 
     return ["using Dune::PDELab::Backend::native;",
             "double maxerror(0.0);",
@@ -905,13 +934,12 @@ def exact_solution_from_dirichlet():
             "  if (std::abs(native({})[i]-native({})[i]) > maxerror)".format(v, solution),
             "    maxerror = std::abs(native({})[i]-native({})[i]);".format(v, solution),
             "",
-            "std::cout << \"Maxerror: \" << maxerror << std::endl;",
+            "std::cout << std::endl << \"Maxerror: \" << maxerror << std::endl;",
             "",
-            "if (maxerror>1e-14)",
+            "if (maxerror>{})".format(get_option("exact_solution_expression")[1]),
             "  throw;"]
 
 
-
 @preamble
 def print_residual():
     ini = name_initree()
@@ -991,8 +1019,9 @@ def vtkoutput():
     dune_solve()
     print_residual()
     print_matrix()
-    if get_option("dirichlet_is_exact_solution"):
-        exact_solution_from_dirichlet()
+
+    if get_option("exact_solution_expression"):
+        exact_solution()
 
     return ["Dune::PDELab::addSolutionToVTKWriter({}, {}, {}, Dune::PDELab::vtk::defaultNameScheme(), {});".format(vtkwriter, gfs, vec, predicate),
             "{}.write({}, Dune::VTK::ascii);".format(vtkwriter, vtkfile)]

test/nonlinear/CMakeLists.txt

@@ -9,6 +9,15 @@ dune_add_system_test(TARGET nonlinear_symdiff
                      SCRIPT dune_vtkcompare.py
                      )
 
+add_generated_executable(UFLFILE nonlinear_dg.ufl
+                         TARGET nonlinear_dg
+                         FORM_COMPILER_ARGS --exact-solution-expression g 1e-2
+                         )
+
+dune_add_system_test(TARGET nonlinear_dg
+                     INIFILE nonlinear_dg.mini
+                     )
+
 # Add the reference code with the PDELab localoperator that produced
 # the reference vtk file
 add_executable(nonlinear_ref reference_main.cc)

test/nonlinear/nonlinear_dg.ufl

+f = Expression("return -2.0*x.size();")
+g = Expression("return x.two_norm2();", on_intersection=True)
+
+cell = triangle
+V = FiniteElement("DG", cell, 1)
+
+u = TrialFunction(V)
+v = TestFunction(V)
+
+n = FacetNormal(cell)('+')
+
+gamma = 1.0
+theta = 1.0
+
+r = inner(grad(u), grad(v))*dx \
+  + inner(n, avg(grad(u)))*jump(v)*dS \
+  + gamma*jump(u)*jump(v)*dS \
+  + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+  - inner(n, grad(u))*v*ds \
+  + gamma*u*v*ds \
+  - theta*u*inner(grad(v), n)*ds \
+  - f*v*dx \
+  + theta*g*inner(grad(v), n)*ds \
+  - gamma*g*v*ds
+
+forms = [r]