From 63006c665b7aeda98f35238ede953899344eeed6 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Wed, 28 Sep 2016 20:26:10 +0200
Subject: [PATCH] First steps towards UFL pullback + implementation of missing
geometry types
---
python/dune/perftool/compile.py | 10 +++++-
python/dune/perftool/loopy/transformer.py | 36 +++++++++++--------
python/dune/perftool/options.py | 14 +++++---
python/dune/perftool/pdelab/argument.py | 10 ------
python/dune/perftool/pdelab/basis.py | 4 +--
python/dune/perftool/pdelab/driver.py | 4 +--
python/dune/perftool/pdelab/geometry.py | 36 ++++++++++++++++++-
python/dune/perftool/pdelab/localoperator.py | 15 ++++++--
python/dune/perftool/pdelab/quadrature.py | 21 -----------
python/dune/perftool/pymbolic/uflmapper.py | 9 ++++-
.../dune/perftool/ufl/modified_terminals.py | 13 +++++++
11 files changed, 111 insertions(+), 61 deletions(-)
diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index 499dcdeb..dfcaeba8 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 bf8fa896..05008cba 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 ef6db4e8..9e345954 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 094d37b8..611f996d 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 23f0a430..912de6ab 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 5bf9b95c..d0b40363 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 8a27c2e1..143a1356 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 9bd174b0..1b65c0a3 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 464251fe..5db88540 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 5eae05a7..c8f10c79 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 82bb8eee..97e4c2c7 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:
--
GitLab