From fe1d653449535ff4d84438a6a3948b189588bb14 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Tue, 6 Sep 2016 17:40:42 +0200
Subject: [PATCH] Make l2 error comparison for systems possible
---
python/dune/perftool/pdelab/driver.py | 170 ++++++++++++++++++++------
python/dune/perftool/ufl/execution.py | 10 +-
test/stokes/CMakeLists.txt | 3 -
test/stokes/stokes.mini | 2 +
test/stokes/stokes_dg.mini | 2 +
test/stokes/stokes_stress.mini | 2 +
test/stokes/stokes_stress.ufl | 7 +-
7 files changed, 154 insertions(+), 42 deletions(-)
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index dfbaf7c4..917ea160 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -807,6 +807,7 @@ def _name_boundary_function(element, boundary, name=None):
define_boundary_function(boundary, name)
return name
+
def name_boundary_function(expr):
# This is a scalar leaf of the ansatz function tree
element, (_, boundary) = get_constraints_metadata(expr)
@@ -814,17 +815,23 @@ def name_boundary_function(expr):
return _name_boundary_function(element, boundary)
-
@symbol
-def name_solution_function(expr):
- # Get expression representing solution
+def name_solution_function(tree_path=()):
+ from dune.perftool.generation import get_global_context_value
+ expr = get_global_context_value("data").object_by_name[get_option("exact_solution_expression")]
+
+ from dune.perftool.ufl.execution import split_expression
+ for i in tree_path:
+ expr = split_expression(expr)[int(i)]
+
from dune.perftool.generation import get_global_context_value
- solution_expression_name = get_option("exact_solution_expression")
- solution_expression = get_global_context_value("data").object_by_name[solution_expression_name]
+ try:
+ name = get_global_context_value("data").object_names[id(expr)]
+ except KeyError:
+ from dune.perftool.generation.cache import _GlobalCounter
+ name = 'generic_{}'.format(_GlobalCounter().get())
- # Create function just like it is done for boundary_function
- name = "solution_function"
- define_boundary_function(solution_expression, name)
+ define_boundary_function(expr, name)
return name
@@ -846,7 +853,7 @@ def interpolate_solution_expression(name):
formdata = _driver_data['formdata']
define_vector(name, formdata)
element = _driver_data['form'].coefficients()[0].ufl_element()
- sol = name_solution_function(element)
+ sol = name_solution_function()
gfs = name_gfs(element)
return "Dune::PDELab::interpolate({}, {}, {});".format(sol,
gfs,
@@ -1181,41 +1188,144 @@ def compare_dofs():
"if (maxerror>{})".format(get_option("compare_dofs")),
" throw;"]
+@symbol
+def type_discrete_grid_function(gfs):
+ return "DGF_{}".format(gfs.upper())
+
@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)]
+ dgf_type = type_discrete_grid_function(gfs)
+ return ["typedef Dune::PDELab::DiscreteGridFunction<decltype({}),decltype({})> {};".format(gfs, vector_name, dgf_type),
+ "{} {}({},{});".format(dgf_type, dgf_name, gfs, vector_name)]
@symbol
def name_discrete_grid_function(gfs, vector_name):
- dgf_name = vector_name + "_dgf"
+ dgf_name = gfs + "_dgf"
define_discrete_grid_function(gfs, vector_name, dgf_name)
return dgf_name
+@symbol
+def type_subgfs(element, tree_path):
+ include_file('dune/pdelab/gridfunctionspace/subspace.hh', filetag='driver')
+ gfs = type_gfs(element)
+ return "Dune::PDELab::GridFunctionSubSpace<{}, Dune::TypeTree::TreePath<{}> >".format(gfs, ', '.join(tree_path))
+
+
@preamble
-def compare_L2_squared():
+def define_subgfs(name, element, tree_path):
+ t = type_subgfs(element, tree_path)
+ gfs = name_gfs(element)
+
+ return '{} {}({});'.format(t, name, gfs)
+
+
+@symbol
+def name_subgfs(element, tree_path):
+ gfs = name_gfs(element)
+ if len(tree_path) == 0:
+ return gfs
+ else:
+ name = "subgfs_{}_{}".format(gfs, '_'.join(tree_path))
+ define_subgfs(name, element, tree_path)
+ return name
+
+
+@preamble
+def typedef_difference_squared_adapter(name, tree_path):
element = _driver_data['form'].coefficients()[0].ufl_element()
formdata = _driver_data['formdata']
- v = name_vector(formdata)
- gfs = name_gfs(element)
- vdgf = name_discrete_grid_function(gfs, v)
- solution_function = name_solution_function(element)
+
+ solution_function = name_solution_function(tree_path)
+ gfs = name_subgfs(element, tree_path)
+ vector = name_vector(formdata)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return 'typedef DifferenceSquaredAdapter<decltype({}), decltype({})> {};'.format(solution_function, dgf, name)
+
+
+@symbol
+def type_difference_squared_adapter(tree_path):
+ t = 'DifferenceSquared_{}'.format('_'.join(tree_path))
+ typedef_difference_squared_adapter(t, tree_path)
+ return t
+
+
+@preamble
+def define_difference_squared_adapter(name, tree_path):
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ formdata = _driver_data['formdata']
+
+ t = type_difference_squared_adapter(tree_path)
+ solution_function = name_solution_function(tree_path)
+ gfs = name_subgfs(element, tree_path)
+ vector = name_vector(formdata)
+ dgf = name_discrete_grid_function(gfs, vector)
+
+ return '{} {}({}, {});'.format(t, name, solution_function, dgf)
+
+
+def name_difference_squared_adapter(tree_path):
+ name = 'differencesquared_{}'.format('_'.join(tree_path))
+ define_difference_squared_adapter(name, tree_path)
+ return name
+
+
+@preamble
+def _accumulate_L2_squared(tree_path):
+ dsa = name_difference_squared_adapter(tree_path)
+ t_dsa = type_difference_squared_adapter(tree_path)
+ accum_error = name_accumulated_l2error()
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")),
+ return ["{",
+ " // L2 error squared of difference between numerical",
+ " // solution and the interpolation of exact solution",
+ " // only subgfs with treepath: {}".format(', '.join(tree_path)),
+ " typename {}::Traits::RangeType err(0.0);".format(t_dsa),
+ " Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
+ " {} += err;".format(accum_error),
+ " std::cout << \"Accumlated L2 Error for tree path {}: \" << err << std::endl;".format(tree_path),
+ "}"]
+
+
+def accumulate_L2_squared(tree_path=()):
+ element = _driver_data['form'].coefficients()[0].ufl_element()
+ from ufl.functionview import select_subelement
+ from ufl.classes import MultiIndex, FixedIndex
+ element = select_subelement(element, MultiIndex(tuple(FixedIndex(int(i)) for i in tree_path)))
+
+ from ufl import MixedElement
+ if isinstance(element, MixedElement):
+ for i, subel in enumerate(element.sub_elements()):
+ accumulate_L2_squared(tree_path + (str(i),))
+ else:
+ _accumulate_L2_squared(tree_path)
+
+
+@preamble
+def define_accumulated_l2error(name):
+ return "double {}(0.0);".format(name)
+
+
+@symbol
+def name_accumulated_l2error():
+ name = 'l2error'
+ define_accumulated_l2error(name)
+ return name
+
+
+@preamble
+def compare_L2_squared():
+ accumulate_L2_squared()
+
+ accum_error = name_accumulated_l2error()
+ return ["std::cout << \"l2errorsquared: \" << {} << std::endl;".format(accum_error),
+ "if ({}>{})".format(accum_error, get_option("compare_l2errorsquared")),
" throw;"]
@@ -1302,10 +1412,6 @@ def vtkoutput():
print_matrix()
if get_option("exact_solution_expression"):
- from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
- raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
-
if get_option("compare_dofs"):
compare_dofs()
if get_option("compare_l2errorsquared"):
@@ -1437,10 +1543,6 @@ def solve_instationary():
print_residual()
print_matrix()
if get_option("exact_solution_expression"):
- from ufl import MixedElement, VectorElement, TensorElement
- if isinstance(_driver_data['form'].coefficients()[0].ufl_element(), (MixedElement, VectorElement, TensorElement)):
- raise NotImplementedError("Comparing to exact solution is olny implemented for scalar elements.")
-
if get_option("compare_dofs"):
compare_dofs()
if get_option("compare_l2errorsquared"):
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index 8bce0a7b..a4b8043f 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -107,8 +107,14 @@ class Expression(Coefficient):
def __mul__(self, other):
# Allow the combination of Expressions
if isinstance(other, Expression):
- from ufl import MixedElement
- return Expression(cppcode=self.c_expr + other.c_expr, element=self.element * other.element)
+ from ufl import MixedElement, TensorElement
+ def get_sub(elem):
+ if isinstance(elem, MixedElement) and not isinstance(elem, (VectorElement, TensorElement)):
+ return elem.sub_elements()
+ else:
+ return [elem]
+ newelem = MixedElement(*(get_sub(self.element) + get_sub(other.element)))
+ return Expression(cppcode=self.c_expr + other.c_expr, element=newelem)
else:
return Coefficient.__mul__(self, other)
diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index 8d6fa96d..0d3751b8 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -5,17 +5,14 @@ dune_symlink_to_source_files(FILES hagenpoiseuille_ref.vtu
dune_add_formcompiler_system_test(UFLFILE stokes.ufl
BASENAME stokes
INIFILE stokes.mini
- SCRIPT dune_vtkcompare.py
)
dune_add_formcompiler_system_test(UFLFILE stokes_dg.ufl
BASENAME stokes_dg
INIFILE stokes_dg.mini
- SCRIPT dune_vtkcompare.py
)
dune_add_formcompiler_system_test(UFLFILE stokes_stress.ufl
BASENAME stokes_stress
INIFILE stokes_stress.mini
- SCRIPT dune_vtkcompare.py
)
diff --git a/test/stokes/stokes.mini b/test/stokes/stokes.mini
index 23a79329..6bb834ff 100644
--- a/test/stokes/stokes.mini
+++ b/test/stokes/stokes.mini
@@ -14,3 +14,5 @@ extension = vtu
[formcompiler]
numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
diff --git a/test/stokes/stokes_dg.mini b/test/stokes/stokes_dg.mini
index 6448f4f5..20d92508 100644
--- a/test/stokes/stokes_dg.mini
+++ b/test/stokes/stokes_dg.mini
@@ -16,3 +16,5 @@ zeroThreshold.data_1 = 1e-6
[formcompiler]
numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-5
diff --git a/test/stokes/stokes_stress.mini b/test/stokes/stokes_stress.mini
index 642ef6d2..4481e91d 100644
--- a/test/stokes/stokes_stress.mini
+++ b/test/stokes/stokes_stress.mini
@@ -14,3 +14,5 @@ extension = vtu
[formcompiler]
numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
diff --git a/test/stokes/stokes_stress.ufl b/test/stokes/stokes_stress.ufl
index 71a32399..d5c41e30 100644
--- a/test/stokes/stokes_stress.ufl
+++ b/test/stokes/stokes_stress.ufl
@@ -1,18 +1,19 @@
v_bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"))
g_p = Expression("8*(1.-x[0])")
-g = g_v * g_p
+g_S = Expression(("0.0", "0.0", "-8*x[1] + 4.", "0.0"))
+g = g_v * g_p * g_S
cell = triangle
P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
P1 = FiniteElement("Lagrange", cell, 1)
-P2_stress = TensorElement("Lagrange", cell, 2)
+P2_stress = TensorElement("DG", cell, 1)
TH = MixedElement(P2, P1, P2_stress)
v, q, T = TestFunctions(TH)
u, p, S = TrialFunctions(TH)
-r = (inner(grad(v), S) + inner(S - grad(u), T)- div(v)*p - q*div(u))*dx
+r = (inner(grad(v), S) + inner(grad(u) - S, T) - div(v)*p - q*div(u))*dx
forms = [r]
--
GitLab