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

@@ -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())

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")

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

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