From da1e31f4daabe378489956c787e08f9314d71438 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Mon, 27 Mar 2017 17:24:48 +0200
Subject: [PATCH] Avoid bug from diagonal_jacobian and splitting for Stokes
We avoid the bug by not applying the transformation and casting the
result of jacobianInverseTransposed to a FieldMatrix. This way we do
not exploit the diagonal structure of the Transformation. Should get
fixed in the near future.
---
python/dune/perftool/options.py | 1 +
python/dune/perftool/pdelab/geometry.py | 14 ++++++++++----
python/dune/perftool/pdelab/localoperator.py | 5 +++--
test/stokes/stokes_quadrilateral.mini | 1 +
4 files changed, 15 insertions(+), 6 deletions(-)
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 9f75b3ce..62da51b6 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -48,6 +48,7 @@ def get_form_compiler_arguments():
parser.add_argument("--sumfact", action="store_true", help="Use sumfactorization")
parser.add_argument("--vectorize-quad", action="store_true", help="whether to generate code with explicit vectorization")
parser.add_argument("--vectorize-grads", action="store_true", help="whether to generate code with explicit vectorization")
+ parser.add_argument("--diagonal_jacobian_bug_workaround", action="store_true", help="Do not use diagonal_jacobian transformation and cast result of jacobianInverseTransposed to FieldMatrx. This is a temporary workaround to avoid a nasty splitting bug.")
# Modify the positional argument to not be a list
args = vars(parser.parse_args())
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 9e98245c..1b9669f4 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -316,10 +316,16 @@ def define_constant_jacobian_inveser_transposed(name, restriction):
globalarg(name, dtype=np.float64, shape=(dim, dim), managed=False)
- return 'auto {} = {}.jacobianInverseTransposed({});'.format(name,
- geo,
- pos,
- )
+ if get_option('diagonal_jacobian_bug_workaround'):
+ jit_type = 'Dune::FieldMatrix<double,{},{}>'.format(dim,dim)
+ else:
+ jit_type = 'auto'
+
+ return '{} {} = {}.jacobianInverseTransposed({});'.format(jit_type,
+ name,
+ geo,
+ pos,
+ )
@backend(interface="define_jit", name="default")
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 088fa742..a925c213 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -376,8 +376,9 @@ def visit_integrals(integrals):
# Maybe make the jacobian inverse diagonal!
if get_option('diagonal_transformation_matrix'):
- from dune.perftool.ufl.transformations.axiparallel import diagonal_jacobian
- integrand = diagonal_jacobian(integrand)
+ if not get_option('diagonal_jacobian_bug_workaround'):
+ from dune.perftool.ufl.transformations.axiparallel import diagonal_jacobian
+ integrand = diagonal_jacobian(integrand)
# Gather dimension indices
from dune.perftool.ufl.dimensionindex import dimension_index_mapping
diff --git a/test/stokes/stokes_quadrilateral.mini b/test/stokes/stokes_quadrilateral.mini
index 0c2a459d..a35a389f 100644
--- a/test/stokes/stokes_quadrilateral.mini
+++ b/test/stokes/stokes_quadrilateral.mini
@@ -14,3 +14,4 @@ extension = vtu
numerical_jacobian = 1, 0 | expand num
exact_solution_expression = g
compare_l2errorsquared = 1e-11
+diagonal_jacobian_bug_workaround = 1
--
GitLab