diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 76aaff919dd3d63045f811ebe9896a5910b529d5..2ac5432e61419d95d75e624952cac09cc107bdec 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -22,7 +22,8 @@ from dune.perftool.pdelab.index import (name_index,
from dune.perftool.pdelab.parameter import (cell_parameter_function,
intersection_parameter_function,
)
-from dune.perftool.pdelab.quadrature import (name_quadrature_weight,
+from dune.perftool.pdelab.quadrature import (pymbolic_quadrature_weight,
+ quadrature_inames,
)
from dune.perftool.pdelab.spaces import (lfs_iname,
)
@@ -102,5 +103,9 @@ class PDELabInterface(object):
# Quadrature related generator functions
#
- def name_quadrature_weight(self):
- return name_quadrature_weight()
+ def pymbolic_quadrature_weight(self):
+ return pymbolic_quadrature_weight()
+
+ # TODO Should this be part of interface or not?
+ def quadrature_inames(self):
+ return quadrature_inames()
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index afcb0e956c1597d526251dcf12efdc86def82d8f..5b72c9b95e54e5796b86c324d6b4c257a15afb06 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -169,11 +169,22 @@ def type_residual():
return "R"
-def name_accumulation_variable(restrictions=(Restriction.NONE,)):
+def name_accumulation_variable(restrictions=None):
ft = get_global_context_value("form_type")
+ measure = get_global_context_value("integral_type")
if ft == 'residual' or ft == 'jacobian_apply':
+ if restrictions is None:
+ if measure == "cell":
+ restrictions = (Restriction.NONE,)
+ else:
+ restrictions = (Restriction.OUTSIDE,)
return name_residual(*restrictions)
if ft == 'jacobian':
+ if restrictions is None:
+ if measure == "cell":
+ restrictions = (Restriction.NONE, Restriction.NONE)
+ else:
+ restrictions = (Restriction.OUTSIDE, Restriction.OUTSIDE)
return name_jacobian(*restrictions)
assert False
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index ad587091096fff8a51debf2002332da2412890bf..aacbcac00980c871c9af99fb73cf537065eb5527 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -1,9 +1,9 @@
""" Generators for basis evaluations """
-from dune.perftool.generation import (backend,
- cached,
+from dune.perftool.generation import (cached,
class_member,
generator_factory,
+ get_backend,
include_file,
instruction,
preamble,
@@ -66,9 +66,8 @@ def evaluate_basis(leaf_element, name, restriction):
temporary_variable(name, shape=(name_lfs_bound(lfs),), decl_method=None)
declare_cache_temporary(leaf_element, restriction, name, which='Function')
cache = name_localbasis_cache(leaf_element)
- qp = name_quadrature_position_in_cell(restriction)
- instruction(inames=(quadrature_iname(),
- ),
+ qp = get_backend("qp_in_cell")(restriction)
+ instruction(inames=get_backend("quad_inames")(),
code='{} = {}.evaluateFunction({}, {}.finiteElement().localBasis());'.format(name,
cache,
qp,
@@ -95,8 +94,7 @@ def evaluate_reference_gradient(leaf_element, name, restriction):
declare_cache_temporary(leaf_element, restriction, name, which='Jacobian')
cache = name_localbasis_cache(leaf_element)
qp = name_quadrature_position_in_cell(restriction)
- instruction(inames=(quadrature_iname(),
- ),
+ instruction(inames=get_backend("quad_inames")(),
code='{} = {}.evaluateJacobian({}, {}.finiteElement().localBasis());'.format(name,
cache,
qp,
@@ -124,7 +122,6 @@ def shape_as_pymbolic(shape):
return tuple(_shape_as_pymbolic(s) for s in shape)
-@backend(interface="eval_coefficient")
@cached
def evaluate_coefficient(element, name, container, restriction, component):
from ufl.functionview import select_subelement
@@ -156,7 +153,7 @@ def evaluate_coefficient(element, name, container, restriction, component):
reduction_expr = Product((coeff, Subscript(Variable(basis), Variable(index))))
instruction(expression=Reduction("sum", index, reduction_expr, allow_simultaneous=True),
assignee=assignee,
- forced_iname_deps=frozenset({quadrature_iname()}).union(frozenset(idims)),
+ forced_iname_deps=frozenset(get_backend("quad_inames")()).union(frozenset(idims)),
forced_iname_deps_is_final=True,
)
@@ -195,6 +192,6 @@ def evaluate_coefficient_gradient(element, name, container, restriction, compone
instruction(expression=Reduction("sum", index, reduction_expr, allow_simultaneous=True),
assignee=assignee,
- forced_iname_deps=frozenset({quadrature_iname()}).union(frozenset(idims)),
+ forced_iname_deps=frozenset(get_backend("quad_inames")()).union(frozenset(idims)),
forced_iname_deps_is_final=True,
)
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 8dd41c1237c43e76d11a10d4b94365d98fbdd472..1210a6a07622a80b90aa2716d9e27515ea8f1148 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -2,6 +2,7 @@ from dune.perftool.ufl.modified_terminals import Restriction
from dune.perftool.pdelab.restriction import restricted_name
from dune.perftool.generation import (cached,
domain,
+ get_backend,
get_global_context_value,
iname,
preamble,
@@ -254,7 +255,7 @@ def define_jacobian_inverse_transposed(name, restriction):
dim = name_dimension()
temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
geo = name_cell_geometry(restriction)
- pos = name_quadrature_position_in_cell(restriction)
+ pos = get_backend("qp_in_cell")(restriction)
return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
geo,
pos,
@@ -273,7 +274,7 @@ def name_jacobian_inverse_transposed(restriction):
def define_jacobian_determinant(name):
temporary_variable(name, shape=())
geo = name_geometry()
- pos = name_quadrature_position()
+ pos = get_backend("quad_pos")()
code = "{} = {}.integrationElement({});".format(name,
geo,
pos,
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 5f9ac0d9bd952d9e21cbba86a4db53477be97944..33af5b7a500dc9627cc9735e306b7baa93513934 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -403,13 +403,13 @@ def generate_accumulation_instruction(visitor, accterm, measure, subdomain_id):
(ansatz_lfs.get_args() + test_lfs.get_args() + (pymbolic_expr,))
)
- from dune.perftool.generation import get_backend, instruction
+ from dune.perftool.generation import instruction
from dune.perftool.options import option_switch
- quad_inames = get_backend(interface="quadinames", selector=option_switch("sumfac"))()
+ quad_inames = visitor.interface.quadrature_inames()
instruction(assignees=(),
expression=expr,
- forced_iname_deps=additional_inames.union(frozenset(visitor.inames).union(quad_inames)),
+ forced_iname_deps=additional_inames.union(frozenset(visitor.inames).union(frozenset(quad_inames))),
forced_iname_deps_is_final=True,
predicates=predicates
)
@@ -459,8 +459,12 @@ def generate_kernel(integrals):
indexmap[j] = indexmap[i]
# Get a transformer instance for this kernel
- from dune.perftool.pdelab import PDELabInterface
- interface = PDELabInterface()
+ if get_option('sumfact'):
+ from dune.perftool.sumfact import SumFactInterface
+ interface = SumFactInterface()
+ else:
+ from dune.perftool.pdelab import PDELabInterface
+ interface = PDELabInterface()
from dune.perftool.ufl.visitor import UFL2LoopyVisitor
visitor = UFL2LoopyVisitor(interface, measure, indexmap)
generate_accumulation_instruction(visitor, term, measure, subdomain_id)
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index b50c8248512ce363cb35b735c199698102ddf110..030c467cc94c14d22bc31c203f9bd2da81953692 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -1,6 +1,7 @@
from dune.perftool.generation import (backend,
cached,
domain,
+ get_backend,
get_global_context_value,
iname,
include_file,
@@ -10,6 +11,8 @@ from dune.perftool.generation import (backend,
from dune.perftool.options import get_option
from dune.perftool.ufl.modified_terminals import Restriction
+from pymbolic.primitives import Variable
+
@iname
def quadrature_iname():
@@ -17,13 +20,13 @@ def quadrature_iname():
return "q"
-@backend(interface="quadinames")
+@backend(interface="quad_inames")
def quadrature_inames():
- return frozenset({quadrature_iname()})
+ return (quadrature_iname(),)
def quadrature_preamble(code, **kw):
- return instruction(inames=quadrature_iname(), code=code, **kw)
+ return instruction(inames=get_backend(interface="quad_inames")(), code=code, **kw)
def name_quadrature_point():
@@ -37,6 +40,7 @@ def define_quadrature_position(name):
)
+@backend(interface="quad_pos")
def name_quadrature_position():
name = "pos"
# To determine the shape, I do query global information here for lack of good alternatives
@@ -52,6 +56,7 @@ def name_quadrature_position():
return name
+@backend(interface="qp_in_cell")
def name_quadrature_position_in_cell(restriction):
if restriction == Restriction.NONE:
return name_quadrature_position()
@@ -67,11 +72,11 @@ def define_quadrature_weight(name):
)
-def name_quadrature_weight():
+def pymbolic_quadrature_weight():
name = 'weight'
temporary_variable(name, shape=())
define_quadrature_weight(name)
- return name
+ return Variable(name)
def estimate_quadrature_order():
diff --git a/python/dune/perftool/pdelab/spaces.py b/python/dune/perftool/pdelab/spaces.py
index 07403a8641ba8701ccd3a89a8c9c619ef8167df9..0e7b0ce29818b72082a55ba795844d2ebf18747f 100644
--- a/python/dune/perftool/pdelab/spaces.py
+++ b/python/dune/perftool/pdelab/spaces.py
@@ -212,6 +212,24 @@ def lfs_iname(element, restriction, count=None, context=''):
return _lfs_iname(element, restriction, context)
+class LFSLocalIndex(FunctionIdentifier):
+ def __init__(self, lfs):
+ self.lfs = lfs
+
+ def __getinitargs__(self):
+ return (self.lfs,)
+
+ @property
+ def name(self):
+ return '{}.local_index'.format(self.lfs)
+
+
+@function_mangler
+def lfs_localindex_mangler(target, func, dtypes):
+ if isinstance(func, LFSLocalIndex):
+ return CallMangleInfo(func.name, (NumpyType(numpy.int32),), (NumpyType(numpy.int32),))
+
+
def name_testfunctionspace(restriction):
return restricted_name("lfsv", restriction)
diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py
index 3f081643c4230fa3168db58beef408640e2a5290..e263c742444e51ccbb301eade94ffcd60472756d 100644
--- a/python/dune/perftool/sumfact/__init__.py
+++ b/python/dune/perftool/sumfact/__init__.py
@@ -1,6 +1,18 @@
-# Trigger some imports that are needed to have all backend implementations visible
-# to the selection mechanisms
-import dune.perftool.sumfact.amatrix
-import dune.perftool.sumfact.sumfact
+from dune.perftool.sumfact.quadrature import (quadrature_inames,
+ quadrature_weight,
+ )
-from dune.perftool.sumfact.sumfact import start_sumfactorization
+from dune.perftool.sumfact.sumfact import pymbolic_trialfunction
+
+from dune.perftool.pdelab import PDELabInterface
+
+
+class SumFactInterface(PDELabInterface):
+ def pymbolic_trialfunction(self, element, restriction, component):
+ return pymbolic_trialfunction(element, restriction, component)
+
+ def quadrature_inames(self):
+ return quadrature_inames()
+
+ def pymbolic_quadrature_weight(self):
+ return quadrature_weight()
diff --git a/python/dune/perftool/sumfact/amatrix.py b/python/dune/perftool/sumfact/amatrix.py
index 968f4f16de04ce7daace2e19f74185ff24380222..1d878add55b3d52c0c539890380e06e3b784a922 100644
--- a/python/dune/perftool/sumfact/amatrix.py
+++ b/python/dune/perftool/sumfact/amatrix.py
@@ -125,6 +125,7 @@ def define_oned_quadrature_weights(name):
def name_oned_quadrature_weights():
name = "qw"
+ globalarg(name, shape=(quadrature_points_per_direction(),), dtype=NumpyType(numpy.float64))
define_oned_quadrature_weights(name)
return name
@@ -138,6 +139,7 @@ def define_oned_quadrature_points(name):
def name_oned_quadrature_points():
name = "qp"
+ globalarg(name, shape=(quadrature_points_per_direction(),), dtype=NumpyType(numpy.float64))
define_oned_quadrature_points(name)
return name
diff --git a/python/dune/perftool/sumfact/quadrature.py b/python/dune/perftool/sumfact/quadrature.py
new file mode 100644
index 0000000000000000000000000000000000000000..ea7104a4c8316f7aff23fd3f4a7f9f6dbdc040e6
--- /dev/null
+++ b/python/dune/perftool/sumfact/quadrature.py
@@ -0,0 +1,114 @@
+from dune.perftool.generation import (backend,
+ domain,
+ function_mangler,
+ get_global_context_value,
+ iname,
+ instruction,
+ temporary_variable,
+ )
+
+from dune.perftool.sumfact.amatrix import (name_number_of_basis_functions_per_direction,
+ name_oned_quadrature_points,
+ name_oned_quadrature_weights,
+ )
+from dune.perftool.pdelab.argument import name_accumulation_variable
+from dune.perftool.pdelab.geometry import dimension_iname
+
+from loopy import CallMangleInfo
+from loopy.symbolic import FunctionIdentifier
+from loopy.types import NumpyType
+
+from pymbolic.primitives import (Call,
+ Product,
+ Subscript,
+ Variable,
+ )
+
+import numpy
+
+
+class BaseWeight(FunctionIdentifier):
+ def __init__(self, accumobj):
+ self.accumobj = accumobj
+
+ def __getinitargs__(self):
+ return (self.accumobj,)
+
+ @property
+ def name(self):
+ return '{}.weight'.format(self.accumobj)
+
+
+@function_mangler
+def base_weight_function_mangler(target, func, dtypes):
+ if isinstance(func, BaseWeight):
+ return CallMangleInfo(func.name, (NumpyType(numpy.float64),), ())
+
+
+@iname
+def sumfact_quad_iname(d, context):
+ name = "quad_{}_{}".format(context, d)
+ domain(name, name_number_of_basis_functions_per_direction())
+ return name
+
+
+@backend(interface="quad_inames", name="sumfact")
+def quadrature_inames(context=''):
+ formdata = get_global_context_value('formdata')
+ dim = formdata.geometric_dimension
+ return tuple(sumfact_quad_iname(d, context) for d in range(dim))
+
+
+def define_recursive_quadrature_weight(name, dir):
+ iname = quadrature_inames()[dir]
+ temporary_variable(name, shape=(), shape_impl=())
+ instruction(expression=Product((recursive_quadrature_weight(dir + 1),
+ Subscript(Variable(name_oned_quadrature_weights()),
+ (Variable(iname),)
+ ))
+ ),
+ assignee=Variable(name),
+ forced_iname_deps=frozenset(quadrature_inames()[dir:]),
+ forced_iname_deps_is_final=True,
+ )
+
+
+def recursive_quadrature_weight(dir=0):
+ formdata = get_global_context_value('formdata')
+ dim = formdata.geometric_dimension
+ if dir == dim:
+ return Call(BaseWeight(name_accumulation_variable()), ())
+ else:
+ name = 'weight_{}'.format(dir)
+ define_recursive_quadrature_weight(name, dir)
+ return Variable(name)
+
+
+def quadrature_weight():
+ return recursive_quadrature_weight()
+
+
+def define_quadrature_position(name):
+ formdata = get_global_context_value('formdata')
+ dim = formdata.geometric_dimension
+ for i in range(dim):
+ instruction(expression=Subscript(Variable(name_oned_quadrature_points()), (Variable(quadrature_inames()[i]),)),
+ assignee=Subscript(Variable(name), (i,)),
+ forced_iname_deps=frozenset(quadrature_inames()),
+ forced_iname_deps_is_final=True,
+ )
+
+
+@backend(interface="quad_pos", name="sumfact")
+def name_quadrature_position():
+ formdata = get_global_context_value('formdata')
+ dim = formdata.geometric_dimension
+ name = 'pos'
+ temporary_variable(name, shape=(dim,), shape_impl=("fv",))
+ define_quadrature_position(name)
+ return name
+
+
+@backend(interface="qp_in_cell", name="sumfact")
+def name_quadrature_position_in_cell(restriction):
+ return name_quadrature_position()
diff --git a/python/dune/perftool/sumfact/sumfact.py b/python/dune/perftool/sumfact/sumfact.py
index 5a3e364d2705281f8184ceb1ab4b013956bba50a..37be6d7d9b7e297181914f54af42a7111608178a 100644
--- a/python/dune/perftool/sumfact/sumfact.py
+++ b/python/dune/perftool/sumfact/sumfact.py
@@ -1,4 +1,6 @@
-from dune.perftool.pdelab.argument import pymbolic_coefficient
+from dune.perftool.pdelab.argument import (name_coefficientcontainer,
+ pymbolic_coefficient,
+ )
from dune.perftool.generation import (backend,
domain,
get_counter,
@@ -23,7 +25,7 @@ from pymbolic.primitives import (Call,
Subscript,
Variable,
)
-
+from dune.perftool.sumfact.quadrature import quadrature_inames
from loopy import Reduction
from pytools import product
@@ -42,7 +44,7 @@ def sumfact_iname(bound, _type):
return name
-def setup_theta(element, container, restriction, component, a_matrices):
+def setup_theta(element, restriction, component, a_matrices):
number_basis = product(mat.cols for mat in a_matrices)
shape = (number_basis,)
inp = get_buffer_temporary("buffer",
@@ -52,6 +54,7 @@ def setup_theta(element, container, restriction, component, a_matrices):
# Write initial coefficients into buffer
basisiname = sumfact_iname(number_basis, "basis")
lfs = name_lfs(element, restriction, component)
+ container = name_coefficientcontainer(restriction)
coeff = pymbolic_coefficient(container, lfs, basisiname)
assignee = Subscript(Variable(inp), (Variable(basisiname),))
return instruction(assignee=assignee,
@@ -60,9 +63,7 @@ def setup_theta(element, container, restriction, component, a_matrices):
)
-# TODO this code is WIP and mainly used for experiments.
-@backend(interface="eval_coefficient", name="sumfact")
-def start_sumfactorization(element, name, container, restriction, component):
+def pymbolic_trialfunction(element, restriction, component):
theta = name_theta()
rows = quadrature_points_per_direction()
cols = basis_functions_per_direction()
@@ -75,16 +76,17 @@ def start_sumfactorization(element, name, container, restriction, component):
num=2
)
- insn_dep = setup_theta(element, container, restriction, component, a_matrices)
+ insn_dep = setup_theta(element, restriction, component, a_matrices)
+ var = sum_factorization_kernel(a_matrices, "buffer", 0, frozenset({insn_dep}))
- sum_factorization_kernel(a_matrices, "buffer", 0, frozenset({insn_dep}))
+ return Subscript(Variable(var), tuple(Variable(i) for i in quadrature_inames()))
- # Do stage 3 (for f=u => mass matrix)
- theta_transposed = name_theta_transposed()
- a_matrix_transposed = AMatrix(theta_transposed, cols, rows)
- a_matrices_transposed = (a_matrix_transposed, a_matrix_transposed)
- return sum_factorization_kernel(a_matrices_transposed, "buffer", 2)
+# # Do stage 3 (for f=u => mass matrix)
+# theta_transposed = name_theta_transposed()
+# a_matrix_transposed = AMatrix(theta_transposed, cols, rows)
+# a_matrices_transposed = (a_matrix_transposed, a_matrix_transposed)
+# var = sum_factorization_kernel(a_matrices_transposed, "buffer", 2)
def sum_factorization_kernel(a_matrices, buffer, stage, insn_dep=frozenset({})):
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index d87c5dcf57f395c6f6ef47359b44825557c776d8..7a7ff15a1eee44714b3d77ac15bad4f04d3a977a 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -256,7 +256,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return Variable(self.interface.name_unit_inner_normal())
def quadrature_weight(self, o):
- return Variable(self.interface.name_quadrature_weight())
+ return self.interface.pymbolic_quadrature_weight()
def jacobian_determinant(self, o):
return Variable(self.interface.name_jacobian_determinant())