diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index 499dcdebde8e4fc6eb20e5ff3a4c4c45639ac337..dfcaeba86bdd21bff5b956aa330f016243e2d7c6 100644
--- a/python/dune/perftool/compile.py
+++ b/python/dune/perftool/compile.py
@@ -36,7 +36,15 @@ def read_ufl(uflfile):
formdatas = []
forms = data.forms
for index, form in enumerate(forms):
- formdatas.append(compute_form_data(form))
+ from dune.perftool.ufl.execution import (JacobianDeterminant,
+ JacobianInverse)
+ formdatas.append(compute_form_data(form,
+ do_apply_function_pullbacks=True,
+ do_apply_geometry_lowering=True,
+ do_apply_integral_scaling=True,
+ preserve_geometry_types=(JacobianDeterminant, JacobianInverse)
+ )
+ )
forms[index] = formdatas[index].preprocessed_form
# We expect at least one form
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index bf8fa8968ce2383558722f7f133985c962e20be6..05008cba2c986a74d761df93dfb9924078be1fbe 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -42,6 +42,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
def __call__(self, o):
# Reset some state variables that are reinitialized for each accumulation term
self.argshape = 0
+ self.transpose_necessary = False
self.redinames = ()
self.inames = []
self.substitution_rules = []
@@ -118,9 +119,6 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
from dune.perftool.pdelab.argument import name_accumulation_variable
accumvar = name_accumulation_variable(arg_restr)
- from dune.perftool.pdelab.quadrature import name_factor
- factor = name_factor()
-
predicates = frozenset({})
# Maybe wrap this instruction into a condiditional. This mostly happens with mixed boundary conditions
@@ -159,9 +157,8 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
predicates = predicates.union(['{} == {}'.format(name, self.subdomain_id)])
from dune.perftool.loopy.functions import PDELabAccumulationFunction
- from pymbolic.primitives import Call, Product
- expr = Product((pymbolic_expr, Variable(factor)))
- expr = Call(PDELabAccumulationFunction(accumvar, len(test_ma)), tuple(a for a in accumargs) + (expr,))
+ from pymbolic.primitives import Call
+ expr = Call(PDELabAccumulationFunction(accumvar, len(test_ma)), tuple(a for a in accumargs) + (pymbolic_expr,))
instruction(assignees=(),
expression=expr,
@@ -199,16 +196,19 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
# For the purpose of basis evaluation, we need to take the leaf element
leaf_element = element.sub_elements()[0]
+ if self.grad:
+ raise ValueError("Gradients should have been transformed to reference gradients!!!")
+
# Have the issued instruction depend on the iname for this localfunction space
iname = lfs_iname(leaf_element, restriction, o.number())
self.inames.append(iname)
- if self.grad:
- from dune.perftool.pdelab.argument import name_testfunction_gradient
- return Subscript(Variable(name_testfunction_gradient(leaf_element, restriction)), (Variable(iname),))
+ if self.reference_grad:
+ from dune.perftool.pdelab.basis import name_reference_gradient
+ return Subscript(Variable(name_reference_gradient(leaf_element, restriction)), (Variable(iname),))
else:
- from dune.perftool.pdelab.argument import name_testfunction
- return Subscript(Variable(name_testfunction(leaf_element, restriction)), (Variable(iname),))
+ from dune.perftool.pdelab.basis import name_basis
+ return Subscript(Variable(name_basis(leaf_element, restriction)), (Variable(iname),))
def coefficient(self, o):
# Do something different for trial function and coefficients from jacobian apply
@@ -219,6 +219,9 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
restriction = Restriction.NEGATIVE
if self.grad:
+ raise ValueError("Gradients should have been transformed to reference gradients!!!")
+
+ if self.reference_grad:
if o.count() == 0:
from dune.perftool.pdelab.argument import name_trialfunction_gradient
return Variable(name_trialfunction_gradient(o.ufl_element(), restriction, self.component))
@@ -262,11 +265,16 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
""" Although indexed is something that we want to handle in the pymbolic mapper we
need to handle it here, as we have to give special treatment to indices into vector
elements: they become a loop index! """
- # Set the generated index for the argument handler to have it available
- self.last_index = self.call(o.ufl_operands[1])
aggr = self.call(o.ufl_operands[0])
+ indices = self.call(o.ufl_operands[1])
+
+ if self.transpose_necessary:
+ assert(len(indices) == 2)
+ assert(self.argshape == 0)
+ indices = (indices[1], indices[0])
+ self.transpose_necessary = False
- use_indices = self.last_index[self.argshape:]
+ use_indices = indices[self.argshape:]
self.argshape = 0
if isinstance(aggr, Subscript):
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index ef6db4e89d56b951b3272397991c237bcb43e6ef..9e345954f24a313d216af9a6061f54edd2e45b81 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -100,8 +100,12 @@ def set_option(key, value):
def get_option(key, default=None):
- """Return the value corresponding to key from option dictionary"""
- # Make sure form compile arguments were read
- init_option_dict()
-
- return _option_dict.get(key, default)
+ try:
+ __IPYTHON__
+ return default
+ except:
+ """Return the value corresponding to key from option dictionary"""
+ # Make sure form compile arguments were read
+ init_option_dict()
+
+ return _option_dict.get(key, default)
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 094d37b8a2f3dad5353b3a168b9b0a84068aa983..611f996d7b21ec8821af5a8ed9aa2b03edcd1fcf 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -17,7 +17,6 @@ from dune.perftool.pdelab.basis import (evaluate_coefficient,
evaluate_coefficient_gradient,
lfs_iname,
name_basis,
- name_basis_gradient,
name_lfs_bound,
)
from dune.perftool import Restriction
@@ -139,15 +138,6 @@ def name_argumentspace(ma):
assert False
-def name_argument(ma):
- if ma.argexpr.number() == 0:
- return name_testfunction(ma)
- if ma.argexpr.number() == 1:
- return name_trialfunction(ma)
- # We should never encounter an argument other than 0 or 1
- assert False
-
-
def name_jacobian(restriction1, restriction2):
# Restrictions may only differ if NONE
if (restriction1 == Restriction.NONE) or (restriction2 == Restriction.NONE):
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 23f0a4300446b42d6206fb202d583faca2d6794e..912de6ab6165b39810a77e59f9e39a822295fef8 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -370,7 +370,6 @@ def evaluate_coefficient(element, name, container, restriction, component):
basis = name_basis(leaf_element, restriction)
if isinstance(sub_element, (VectorElement, TensorElement)):
- # TODO we need a concept of symmetry here
lfs = lfs_child(lfs, idims, shape=shape_as_pymbolic(shape), symmetry=element.symmetry)
from dune.perftool.pdelab.argument import pymbolic_coefficient
@@ -406,10 +405,9 @@ def evaluate_coefficient_gradient(element, name, container, restriction, compone
temporary_variable(name, shape=shape, shape_impl=shape_impl)
lfs = name_lfs(element, restriction, component)
index = lfs_iname(leaf_element, restriction, context='trialgrad')
- basis = name_basis_gradient(leaf_element, restriction)
+ basis = name_reference_gradient(leaf_element, restriction)
if isinstance(sub_element, (VectorElement, TensorElement)):
- # TODO we need a concept of symmetry here
lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]), symmetry=element.symmetry)
from dune.perftool.pdelab.argument import pymbolic_coefficient
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 5bf9b95c949e46abf6e88c9e982d0879c3ecedb2..d0b403633da5c79610b261be51a8242446f44a73 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -1076,7 +1076,7 @@ def name_timing_stream():
@preamble
def dune_solve():
# Test if form is linear in ansatzfunction
- linear = is_linear(_driver_data['form'])
+ linear = is_linear(_driver_data['formdata'].original_form)
# Test wether we want to do matrix free operator evaluation
matrix_free = get_option('matrix_free')
@@ -1499,7 +1499,7 @@ def time_loop():
def solve_instationary():
# Test if form is linear in ansatzfunction
- linear = is_linear(_driver_data['form'])
+ linear = is_linear(_driver_data['formdata'].original_form)
# Test wether we want to do matrix free operator evaluation
matrix_free = get_option('matrix_free')
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 8a27c2e179e5e8a7bc4cb568ce258066612853a3..143a1356214859963345328a0529c07e9800156d 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -11,6 +11,7 @@ from dune.perftool.pdelab.quadrature import (name_quadrature_position,
quadrature_preamble,
)
from ufl.algorithms import MultiFunction
+from pymbolic.primitives import Variable
class GeometryMapper(MultiFunction):
@@ -38,10 +39,23 @@ class GeometryMapper(MultiFunction):
# TODO This one was just copied over so, it might need some tweaking
def facet_area(self, o):
- from pymbolic.primitives import Variable
from dune.perftool.pdelab.geometry import name_facetarea
return Variable(name_facetarea())
+ def quadrature_weight(self, o):
+ from dune.perftool.pdelab.quadrature import name_quadrature_weight
+ return Variable(name_quadrature_weight())
+
+ def jacobian_determinant(self, o):
+ return Variable(name_jacobian_determinant())
+
+ def jacobian_inverse(self, o):
+ self.transpose_necessary = True
+ return Variable(name_jacobian_inverse_transposed(self.restriction))
+
+ def jacobian(self, o):
+ raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
+
@iname
def _dimension_iname(context, count):
@@ -307,3 +321,23 @@ def name_jacobian_inverse_transposed(restriction):
name = restricted_name("jit", restriction)
define_jacobian_inverse_transposed(name, restriction)
return name
+
+
+def define_jacobian_determinant(name):
+ temporary_variable(name, shape=())
+ geo = name_geometry()
+ pos = name_quadrature_position()
+ code = "{} = {}.integrationElement({});".format(name,
+ geo,
+ pos,
+ )
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ read_variables=frozenset({pos}),
+ )
+
+
+def name_jacobian_determinant():
+ name = 'detjac'
+ define_jacobian_determinant(name)
+ return name
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 9bd174b0f7abff14c71bdf2e0159a1e94544be21..1b65c0a33653917aa1065d86789fdbe7e86cb799 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -187,7 +187,7 @@ def assembly_routine_signature(formdata):
# Check if form is linear
from dune.perftool.pdelab.driver import is_linear
- linear = is_linear(formdata.preprocessed_form)
+ linear = is_linear(formdata.original_form)
if form_type == 'residual':
if integral_type == 'cell':
@@ -509,8 +509,17 @@ def generate_localoperator_kernels(formdata, data):
# Generate the necessary jacobian methods
if not get_option("numerical_jacobian"):
from ufl import derivative
- from ufl.algorithms import expand_derivatives
- jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
+ from ufl.algorithms import compute_form_data
+ from dune.perftool.ufl.execution import (JacobianDeterminant,
+ JacobianInverse,
+ )
+ jacform = derivative(formdata.original_form, formdata.original_form.coefficients()[0])
+ jacform = compute_form_data(jacform,
+ do_apply_function_pullbacks=True,
+ do_apply_geometry_lowering=True,
+ do_apply_integral_scaling=True,
+ preserve_geometry_types=(JacobianDeterminant, JacobianInverse),
+ ).preprocessed_form
with global_context(form_type="jacobian"):
for measure in set(i.integral_type() for i in jacform.integrals()):
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 464251fe33bceb7b528cd7cce8eabcff79ad4b42..5db88540ad6e001d6a7753414924e17ebc543b8d 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -57,27 +57,6 @@ def name_quadrature_weight():
return "{}.weight()".format(qp)
-def define_quadrature_factor(fac):
- weight = name_quadrature_weight()
- from dune.perftool.pdelab.geometry import name_geometry
- geo = name_geometry()
- pos = name_quadrature_position()
- code = "{} = {}*{}.integrationElement({});".format(fac,
- weight,
- geo,
- pos,
- )
- return quadrature_preamble(code,
- assignees=frozenset({fac}),
- read_variables=frozenset({pos}),
- )
-
-
-def name_factor():
- temporary_variable("fac", shape=())
- define_quadrature_factor("fac")
- return "fac"
-
def name_order():
# TODO how to determine integration order in a generic fashion here?
diff --git a/python/dune/perftool/pymbolic/uflmapper.py b/python/dune/perftool/pymbolic/uflmapper.py
index 5eae05a7121eed8ddc692f66052cc80868bf55a7..c8f10c797d698ac2ce1e8305c5c04dcd986d0a96 100644
--- a/python/dune/perftool/pymbolic/uflmapper.py
+++ b/python/dune/perftool/pymbolic/uflmapper.py
@@ -8,7 +8,7 @@ This is mainly intended as a base class for anything mapping ufl to pymbolic.
from __future__ import absolute_import
from ufl.algorithms import MultiFunction
-from pymbolic.primitives import Product, Quotient, Subscript, Sum, Variable
+from pymbolic.primitives import Product, Quotient, Subscript, Sum, Variable, Call
# The constructed pymbolic expressions use n-ary operators instead of ufls binary operators
from dune.perftool.ufl.flatoperators import get_operands
@@ -39,3 +39,10 @@ class UFL2PymbolicMapper(MultiFunction):
def zero(self, o):
return 0
+
+ def abs(self, o):
+ from ufl.classes import JacobianDeterminant
+ if isinstance(o.ufl_operands[0], JacobianDeterminant):
+ return self.call(o.ufl_operands[0])
+ else:
+ return Call('abs', self.call(o.ufl_operands[0]))
diff --git a/python/dune/perftool/ufl/modified_terminals.py b/python/dune/perftool/ufl/modified_terminals.py
index 82bb8eee94b2c0ec9e3a46835faae0e0dc610e70..97e4c2c75e1a8171cfbf9bb6f192763e8f983b65 100644
--- a/python/dune/perftool/ufl/modified_terminals.py
+++ b/python/dune/perftool/ufl/modified_terminals.py
@@ -13,6 +13,7 @@ class ModifiedTerminalTracker(MultiFunction):
def __init__(self):
MultiFunction.__init__(self)
self.grad = False
+ self.reference = True
self.reference_grad = False
self.restriction = Restriction.NONE
self.component = MultiIndex(())
@@ -51,12 +52,19 @@ class ModifiedTerminalTracker(MultiFunction):
self.component = MultiIndex(())
return ret
+ def reference_value(self, o):
+ self.reference = True
+ ret = self.call(o.ufl_operands[0])
+ self.reference = False
+ return ret
+
class ModifiedArgumentDescriptor(MultiFunction):
def __init__(self, e):
MultiFunction.__init__(self)
self.grad = False
+ self.reference = False
self.reference_grad = False
self.index = None
self.restriction = Restriction.NONE
@@ -77,6 +85,10 @@ class ModifiedArgumentDescriptor(MultiFunction):
self.reference_grad = True
self(o.ufl_operands[0])
+ def reference_value(self, o):
+ self.reference = True
+ self(o.ufl_operands[0])
+
def positive_restricted(self, o):
self.restriction = Restriction.POSITIVE
self(o.ufl_operands[0])
@@ -129,6 +141,7 @@ class _ModifiedArgumentExtractor(MultiFunction):
grad = pass_on
reference_grad = pass_on
function_view = pass_on
+ reference_value = pass_on
def argument(self, o):
if self.argnumber is None or o.number() == self.argnumber: