From d1c317e55c9cf679916952b5cbf76e5118c0bc71 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Wed, 10 Aug 2016 14:46:03 +0200
Subject: [PATCH] Formcompiler argument marking exact solution
NOTE: Works only for scalar equations
---
python/dune/perftool/options.py | 6 ++--
python/dune/perftool/pdelab/driver.py | 47 ++++++++++++++++++++++-----
test/nonlinear/CMakeLists.txt | 9 +++++
test/nonlinear/nonlinear_dg.ufl | 26 +++++++++++++++
4 files changed, 76 insertions(+), 12 deletions(-)
create mode 100644 test/nonlinear/nonlinear_dg.ufl
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 1aeb7f27..e4b51c3a 100644
--- a/python/dune/perftool/options.py
+++ b/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
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index a20571f8..f7a27ff0 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/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)]
diff --git a/test/nonlinear/CMakeLists.txt b/test/nonlinear/CMakeLists.txt
index ff8e1624..5ebedaa8 100644
--- a/test/nonlinear/CMakeLists.txt
+++ b/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)
diff --git a/test/nonlinear/nonlinear_dg.ufl b/test/nonlinear/nonlinear_dg.ufl
new file mode 100644
index 00000000..bf606f9c
--- /dev/null
+++ b/test/nonlinear/nonlinear_dg.ufl
@@ -0,0 +1,26 @@
+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]
--
GitLab