diff --git a/dune/perftool/localbasiscache.hh b/dune/perftool/localbasiscache.hh
index 8d3449105944560ac6dd8058f59b3d34bba38bd0..74896ec9dbbf17390ecb57200f8a73a1658c981d 100644
--- a/dune/perftool/localbasiscache.hh
+++ b/dune/perftool/localbasiscache.hh
@@ -60,4 +60,28 @@ class LocalBasisCacheWithoutReferences
Cache c;
};
+
+template <typename G, typename T>
+void fillQuadratureWeightsCache(const G& geo, const int quadOrder, T& quadratureWeights){
+ if(quadratureWeights.size() != 0)
+ return;
+ else{
+ for (const auto& qp : quadratureRule(geo, quadOrder)){
+ quadratureWeights.push_back(qp.weight());
+ }
+ }
+}
+
+
+template <typename G, typename T>
+void fillQuadraturePointsCache(const G& geo, const int quadOrder, T& quadraturePoints){
+ if(quadraturePoints.size() != 0)
+ return;
+ else{
+ for (const auto& qp : quadratureRule(geo, quadOrder)){
+ quadraturePoints.push_back(qp.position());
+ }
+ }
+}
+
#endif
diff --git a/python/dune/perftool/loopy/target.py b/python/dune/perftool/loopy/target.py
index 97ba3c01ae8ba5f816e5126de3912cb0a7909462..8854688164e7805d1de90df31aa157b85c0efca6 100644
--- a/python/dune/perftool/loopy/target.py
+++ b/python/dune/perftool/loopy/target.py
@@ -46,18 +46,6 @@ class DuneASTBuilder(CASTBuilder):
if temp_var.decl_method:
return Line(temp_var.decl_method(temp_var.name, temp_var.shape, temp_var.shape_impl))
- def emit_sequential_loop(self, codegen_state, iname, iname_dtype, static_lbound, static_ubound, inner):
- # Some of our loops are special as they use PDELab specific concepts.
- # Fortunately those loops are tied to specific inames.
- from dune.perftool.pdelab.quadrature import quadrature_iname
- if iname == quadrature_iname():
- from cgen import CustomLoop
- from dune.perftool.pdelab.quadrature import quadrature_loop_statement
- loop_stmt = quadrature_loop_statement()
- return CustomLoop("for({})".format(loop_stmt), inner)
- else:
- return CASTBuilder.emit_sequential_loop(self, codegen_state, iname, iname_dtype, static_lbound, static_ubound, inner)
-
class DuneTarget(TargetBase):
def __init__(self):
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index 17d967b1c2ee36e732fe790fba7ebe92f3ec8624..42c475567e348aa57b08198bad54ef8fc5872599 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -16,9 +16,6 @@ from dune.perftool.pdelab.spaces import (lfs_child,
name_lfs_bound,
type_gfs,
)
-from dune.perftool.pdelab.quadrature import (name_quadrature_position_in_cell,
- quadrature_iname,
- )
from dune.perftool.pdelab.geometry import (dimension_iname,
name_dimension,
name_jacobian_inverse_transposed,
@@ -27,6 +24,7 @@ from dune.perftool.pdelab.geometry import (dimension_iname,
from dune.perftool.pdelab.localoperator import (lop_template_ansatz_gfs,
lop_template_test_gfs,
)
+from dune.perftool.tools import get_pymbolic_basename
from dune.perftool.pdelab.driver import FEM_name_mangling
from dune.perftool.pdelab.restriction import restricted_name
from pymbolic.primitives import Product, Subscript, Variable
@@ -71,11 +69,11 @@ def evaluate_basis(leaf_element, name, restriction):
instruction(inames=get_backend("quad_inames")(),
code='{} = {}.evaluateFunction({}, {}.finiteElement().localBasis());'.format(name,
cache,
- qp,
+ str(qp),
lfs,
),
assignees=frozenset({name}),
- read_variables=frozenset({qp}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
)
@@ -93,15 +91,15 @@ def evaluate_reference_gradient(leaf_element, name, restriction):
lfs = name_leaf_lfs(leaf_element, restriction)
temporary_variable(name, shape=(name_lfs_bound(lfs), 1, name_dimension()), decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian'))
cache = name_localbasis_cache(leaf_element)
- qp = name_quadrature_position_in_cell(restriction)
+ qp = get_backend("qp_in_cell")(restriction)
instruction(inames=get_backend("quad_inames")(),
code='{} = {}.evaluateJacobian({}, {}.finiteElement().localBasis());'.format(name,
cache,
- qp,
+ str(qp),
lfs,
),
assignees=frozenset({name}),
- read_variables=frozenset({qp}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
)
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 1210a6a07622a80b90aa2716d9e27515ea8f1148..63ffcd08ad8362c276b54034f9bc4a4e7a74ca47 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -5,15 +5,43 @@ from dune.perftool.generation import (cached,
get_backend,
get_global_context_value,
iname,
+ include_file,
preamble,
temporary_variable,
)
-from dune.perftool.pdelab.quadrature import (name_quadrature_position,
- name_quadrature_position_in_cell,
+from dune.perftool.pdelab.quadrature import (pymbolic_quadrature_position_in_cell,
quadrature_preamble,
)
+from dune.perftool.tools import get_pymbolic_basename
from ufl.algorithms import MultiFunction
from pymbolic.primitives import Variable
+from pymbolic.primitives import Expression as PymbolicExpression
+
+
+@preamble
+def define_reference_element(name):
+ geo = name_geometry()
+ include_file("dune/pdelab/common/referenceelements.hh", filetag="operatorfile")
+ return "auto {} = referenceElement({});".format(name, geo)
+
+
+def name_reference_element():
+ name = "refElement"
+ define_reference_element(name)
+ return name
+
+
+@preamble
+def define_localcenter(name):
+ reference_element = name_reference_element()
+ # Note: position(i,c) stands for the barycenter of the i-th subentity of codimension c"
+ return "auto {} = {}.position(0,0);".format(name, reference_element)
+
+
+def name_localcenter():
+ name = "localcenter"
+ define_localcenter(name)
+ return name
@iname
@@ -167,25 +195,27 @@ def apply_in_cell_transformation(name, local, restriction):
geo = name_in_cell_geometry(restriction)
return quadrature_preamble("{} = {}.global({});".format(name,
geo,
- local,
+ str(local),
),
assignees=frozenset({name}),
- read_variables=frozenset({local}),
+ read_variables=frozenset({get_pymbolic_basename(local)}),
)
-def name_in_cell_coordinates(local, basename, restriction):
+def pymbolic_in_cell_coordinates(local, restriction):
+ basename = get_pymbolic_basename(local)
name = "{}_in_{}side".format(basename, "in" if restriction is Restriction.NEGATIVE else "out")
temporary_variable(name, shape=(name_dimension(),), shape_impl=("fv",))
apply_in_cell_transformation(name, local, restriction)
- return name
+ return Variable(name)
-def to_cell_coordinates(local, basename, restriction):
+def to_cell_coordinates(local, restriction):
+ assert isinstance(local, PymbolicExpression)
if restriction == Restriction.NONE:
return local
else:
- return name_in_cell_coordinates(local, basename, restriction)
+ return pymbolic_in_cell_coordinates(local, restriction)
def name_dimension():
@@ -200,10 +230,10 @@ def name_intersection_dimension():
def evaluate_unit_outer_normal(name):
ig = name_intersection_geometry_wrapper()
- qp = name_quadrature_position()
+ qp = get_backend("quad_pos")()
return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
assignees=frozenset({name}),
- read_variables=frozenset({qp}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
)
@@ -258,10 +288,10 @@ def define_jacobian_inverse_transposed(name, restriction):
pos = get_backend("qp_in_cell")(restriction)
return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
geo,
- pos,
+ str(pos),
),
assignees=frozenset({name}),
- read_variables=frozenset({pos}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
)
@@ -277,11 +307,11 @@ def define_jacobian_determinant(name):
pos = get_backend("quad_pos")()
code = "{} = {}.integrationElement({});".format(name,
geo,
- pos,
+ str(pos),
)
return quadrature_preamble(code,
assignees=frozenset({name}),
- read_variables=frozenset({pos}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
)
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 0d0b06d1b9aa670457cbe1a49bc62ae8092e0eb5..48998b88d7ba14d14e70f62cba9436f0d670eeb1 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -508,12 +508,6 @@ def generate_kernel(integrals):
from dune.perftool.loopy.stages import finalize_stage_instructions
kernel = finalize_stage_instructions(kernel)
- # This is also silly, but 1D DG Schemes never need the geometry, so the quadrature
- # statement actually introduces a preamble at a stage where preambles cannot be generated
- # anymore. TODO think about how to avoid this!!!
- from dune.perftool.pdelab.quadrature import quadrature_loop_statement
- quadrature_loop_statement()
-
# Now add the preambles to the kernel
preambles = [(i, p) for i, p in enumerate(retrieve_cache_items("preamble"))]
kernel = kernel.copy(preambles=preambles)
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index cdb57bdea95ad4c6b7719b59826fea6311ebe512..602182163e0ad531e7992a04525a9cc36acd4802 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -5,6 +5,7 @@ from dune.perftool.generation import (cached,
class_member,
constructor_parameter,
generator_factory,
+ get_backend,
initializer_list,
preamble,
temporary_variable
@@ -13,10 +14,11 @@ from dune.perftool.pdelab.geometry import (name_cell,
name_dimension,
name_intersection,
)
-from dune.perftool.pdelab.quadrature import (name_quadrature_position,
- name_quadrature_position_in_cell,
+from dune.perftool.pdelab.quadrature import (pymbolic_quadrature_position,
+ pymbolic_quadrature_position_in_cell,
quadrature_preamble,
)
+from dune.perftool.tools import get_pymbolic_basename
from dune.perftool.cgen.clazz import AccessModifier
from dune.perftool.pdelab.localoperator import (class_type_from_cache,
localoperator_basename,
@@ -123,20 +125,21 @@ def define_parameter_function_class_member(name, expr, baset, shape, cell):
@preamble
def evaluate_cellwise_constant_parameter_function(name, restriction):
- dim = name_dimension()
param = name_paramclass()
entity = name_cell(restriction)
+ from dune.perftool.pdelab.geometry import name_localcenter
+ pos = name_localcenter()
from dune.perftool.generation.loopy import valuearg
import numpy
valuearg(name, dtype=numpy.float64)
- return 'auto {} = {}.{}({}, Dune::FieldVector<double, {}>(0.0));'.format(name,
- name_paramclass(),
- name,
- entity,
- dim,
- )
+ return 'auto {} = {}.{}({}, {});'.format(name,
+ name_paramclass(),
+ name,
+ entity,
+ pos,
+ )
@preamble
@@ -145,35 +148,35 @@ def evaluate_intersectionwise_constant_parameter_function(name):
from dune.perftool.generation import get_global_context_value
it = get_global_context_value("integral_type")
assert it is not 'cell'
- dim = name_dimension()
param = name_paramclass()
intersection = name_intersection()
+ pos = name_localcenter()
from dune.perftool.generation.loopy import valuearg
import numpy
valuearg(name, dtype=numpy.float64)
- return 'auto {} = {}.{}({}, Dune::FieldVector<double, {}>(0.0));'.format(name,
- name_paramclass(),
- name,
- intersection,
- dim,
- )
+ return 'auto {} = {}.{}({}, {});'.format(name,
+ name_paramclass(),
+ name,
+ intersection,
+ pos,
+ )
def evaluate_cell_parameter_function(name, restriction):
param = name_paramclass()
entity = name_cell(restriction)
- pos = name_quadrature_position_in_cell(restriction)
+ pos = get_backend(interface="qp_in_cell")(restriction)
return quadrature_preamble('{} = {}.{}({}, {});'.format(name,
name_paramclass(),
name,
entity,
- pos,
+ str(pos),
),
assignees=frozenset({name}),
- read_variables=frozenset({pos}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
)
@@ -185,15 +188,15 @@ def evaluate_intersection_parameter_function(name):
param = name_paramclass()
intersection = name_intersection()
- pos = name_quadrature_position()
+ pos = get_backend("quad_pos")()
return quadrature_preamble('{} = {}.{}({}, {});'.format(name,
name_paramclass(),
name,
intersection,
- pos,
+ str(pos),
),
assignees=frozenset({name}),
- read_variables=frozenset({pos}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
)
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 030c467cc94c14d22bc31c203f9bd2da81953692..f00e2e8f3062163a7a24fb2d7fac3ec97b62bd0a 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -1,22 +1,60 @@
+import numpy
+
from dune.perftool.generation import (backend,
cached,
+ class_member,
domain,
get_backend,
get_global_context_value,
+ globalarg,
iname,
include_file,
instruction,
+ preamble,
temporary_variable,
+ valuearg,
)
+from dune.perftool.pdelab.localoperator import lop_template_range_field
from dune.perftool.options import get_option
from dune.perftool.ufl.modified_terminals import Restriction
-from pymbolic.primitives import Variable
+from pymbolic.primitives import Variable, Subscript
+
+
+@preamble
+def define_quadrature_rule(name):
+ include_file('dune/pdelab/common/quadraturerules.hh', filetag='operatorfile')
+ from dune.perftool.pdelab.geometry import name_geometry
+ geo = name_geometry()
+ order = name_quadrature_order()
+ return "const auto {} = quadratureRule({}, {});".format(name, geo, order)
+
+
+def name_quadrature_rule():
+ name = "quadrature_rule"
+ define_quadrature_rule(name)
+ return name
+
+
+@preamble
+def define_quadrature_bound(name):
+ quad_rule = name_quadrature_rule()
+ return "auto {} = {}.size();".format(name, quad_rule)
+
+
+def name_quadrature_bound():
+ name = "quadrature_size"
+ define_quadrature_bound(name)
+
+ # Quadrature bound is a valuearg for loopy
+ valuearg(name, dtype=numpy.int32)
+
+ return name
@iname
def quadrature_iname():
- domain("q", "q_n")
+ domain("q", name_quadrature_bound())
return "q"
@@ -30,56 +68,123 @@ def quadrature_preamble(code, **kw):
def name_quadrature_point():
+ # Note: Used for qp_in_inside/qp_in_outside
return "qp"
-def define_quadrature_position(name):
- qp = name_quadrature_point()
- return quadrature_preamble(code="{} = {}.position();".format(name, qp),
- assignees=frozenset({name}),
- )
-
-
-@backend(interface="quad_pos")
-def name_quadrature_position():
- name = "pos"
+def _local_dim():
# To determine the shape, I do query global information here for lack of good alternatives
from dune.perftool.generation import get_global_context_value
it = get_global_context_value("integral_type")
from dune.perftool.pdelab.geometry import name_dimension, name_intersection_dimension
if it == 'cell':
- shape = (name_dimension(),)
+ dim = name_dimension()
else:
- shape = (name_intersection_dimension(),)
- temporary_variable(name, shape=shape, shape_impl=("fv",))
- define_quadrature_position(name)
+ dim = name_intersection_dimension()
+
+ return dim
+
+
+@preamble
+def fill_quadrature_points_cache(name):
+ from dune.perftool.pdelab.geometry import name_geometry
+ geo = name_geometry()
+ quad_order = quadrature_order()
+ include_file("dune/perftool/localbasiscache.hh", filetag='operatorfile')
+ return "fillQuadraturePointsCache({}, {}, {});".format(geo, quad_order, name)
+
+
+@class_member("operator")
+def typedef_quadrature_points(name):
+ range_field = lop_template_range_field()
+ dim = _local_dim()
+ return "using {} = typename Dune::QuadraturePoint<{}, {}>::Vector;".format(name, range_field, dim)
+
+
+def type_quadrature_points(name):
+ name = name.upper()
+ typedef_quadrature_points(name)
return name
+@class_member("operator")
+def define_quadrature_points(name):
+ qp_type = type_quadrature_points(name)
+ return "mutable std::vector<{}> {};".format(qp_type, name)
+
+
+def name_quadrature_points():
+ """Name of vector storing quadrature points as class member"""
+ dim = _local_dim()
+ name = "qp_order" + str(dim)
+ define_quadrature_points(name)
+ fill_quadrature_points_cache(name)
+ return name
+
+
+@backend(interface="quad_pos")
+def pymbolic_quadrature_position():
+ quad_points = name_quadrature_points()
+ quad_iname = quadrature_iname()
+ from pymbolic.primitives import Subscript, Variable
+ return Subscript(Variable(quad_points), (Variable(quad_iname),))
+
+
@backend(interface="qp_in_cell")
-def name_quadrature_position_in_cell(restriction):
- if restriction == Restriction.NONE:
- return name_quadrature_position()
- else:
- from dune.perftool.pdelab.geometry import to_cell_coordinates
- return to_cell_coordinates(name_quadrature_position(), name_quadrature_point(), restriction)
+def pymbolic_quadrature_position_in_cell(restriction):
+ from dune.perftool.pdelab.geometry import to_cell_coordinates
+ return to_cell_coordinates(pymbolic_quadrature_position(), restriction)
+
+@preamble
+def fill_quadrature_weights_cache(name):
+ from dune.perftool.pdelab.geometry import name_geometry
+ geo = name_geometry()
+ quad_order = quadrature_order()
+ include_file("dune/perftool/localbasiscache.hh", filetag='operatorfile')
+ return "fillQuadratureWeightsCache({}, {}, {});".format(geo, quad_order, name)
-def define_quadrature_weight(name):
- qp = name_quadrature_point()
- return quadrature_preamble(code="auto {} = {}.weight();".format(name, qp),
- assignees=frozenset({name})
- )
+
+@class_member("operator")
+def typedef_quadrature_weights(name):
+ range_field = lop_template_range_field()
+ dim = _local_dim()
+ return "using {} = typename Dune::QuadraturePoint<{}, {}>::Field;".format(name, range_field, dim)
def pymbolic_quadrature_weight():
- name = 'weight'
- temporary_variable(name, shape=())
- define_quadrature_weight(name)
- return Variable(name)
+ vec = name_quadrature_weights()
+ return Subscript(Variable(vec),
+ tuple(Variable(i) for i in quadrature_inames()))
+
+
+def type_quadrature_weights(name):
+ name = name.upper()
+ typedef_quadrature_weights(name)
+ return name
+
+
+@class_member("operator")
+def define_quadrature_weights(name):
+ qw_type = type_quadrature_weights(name)
+ return "mutable std::vector<{}> {};".format(qw_type, name)
+
+
+def name_quadrature_weights():
+ """"Name of vector storing quadrature weights as class member"""
+ dim = _local_dim()
+ name = "qw_order" + str(dim)
+ define_quadrature_weights(name)
+ fill_quadrature_weights_cache(name)
+
+ # Quadrature weighs is a globar argument for loopy
+ shape = name_quadrature_bound()
+ globalarg(name, shape=(shape,), dtype=numpy.float64)
+
+ return name
-def estimate_quadrature_order():
+def _estimate_quadrature_order():
"""Estimate quadrature order using polynomial degree estimation from UFL"""
# According to UFL documentation estimate_total_polynomial_degree
# should only be called on preprocessed forms.
@@ -92,32 +197,24 @@ def estimate_quadrature_order():
for i in integrals:
polynomial_degree = max(polynomial_degree, i.metadata()['estimated_polynomial_degree'])
- # Note: In PDELab quadrature order m means that integration of
- # polynomials of degree m is exact.
return polynomial_degree
-def name_order():
- # Quadrature order from UFL estimation
- quadrature_order = estimate_quadrature_order()
+def quadrature_order():
+ """Get quadrature order
- # If set use quadrature_order from formcompiler arguments
- #
- # TODO: make it possible to choose different quadrature order for
- # different integral types?
+ Note: In PDELab quadrature order m means that integration of
+ polynomials of degree m is excat.
+ """
if get_option("quadrature_order"):
quadrature_order = get_option("quadrature_order")
+ else:
+ quadrature_order = _estimate_quadrature_order()
- return str(quadrature_order)
+ return quadrature_order
-def quadrature_loop_statement():
- include_file('dune/pdelab/common/quadraturerules.hh', filetag='operatorfile')
- qp = name_quadrature_point()
- from dune.perftool.pdelab.geometry import name_geometry
- geo = name_geometry()
- order = name_order()
- return "const auto& {} : quadratureRule({}, {})".format(qp,
- geo,
- order,
- )
+def name_quadrature_order():
+ # Quadrature order from UFL estimation
+ quad_order = quadrature_order()
+ return str(quad_order)
diff --git a/python/dune/perftool/sumfact/amatrix.py b/python/dune/perftool/sumfact/amatrix.py
index 1d878add55b3d52c0c539890380e06e3b784a922..8dfc180cdb10d6efac8d6089335487ea847f25b7 100644
--- a/python/dune/perftool/sumfact/amatrix.py
+++ b/python/dune/perftool/sumfact/amatrix.py
@@ -20,7 +20,7 @@ from dune.perftool.generation import (class_member,
from dune.perftool.pdelab.localoperator import (name_domain_field,
lop_template_range_field,
)
-from dune.perftool.pdelab.quadrature import estimate_quadrature_order
+from dune.perftool.pdelab.quadrature import quadrature_order
from loopy import CallMangleInfo
from loopy.symbolic import FunctionIdentifier
@@ -75,7 +75,7 @@ def quadrature_points_per_direction():
# as soon as we have a generic implementation
# Quadrature order
- q = estimate_quadrature_order()
+ q = quadrature_order()
# Quadrature points in per direction
nb_qp = q // 2 + 1
@@ -232,7 +232,7 @@ def define_theta(name, transpose):
shape = (basis_functions_per_direction(), quadrature_points_per_direction())
else:
shape = (quadrature_points_per_direction(), basis_functions_per_direction())
- globalarg(name, shape=shape, dtype=numpy.float32, dim_tags="f,f")
+ globalarg(name, shape=shape, dtype=numpy.float64, dim_tags="f,f")
initializer_list(name, [str(axis) for axis in shape], classtag="operator")
construct_theta(name, transpose)
return "{} {};".format(theta_type, name)
diff --git a/python/dune/perftool/tools.py b/python/dune/perftool/tools.py
new file mode 100644
index 0000000000000000000000000000000000000000..f528afcd8459538d49f7777429fa52dd6318b253
--- /dev/null
+++ b/python/dune/perftool/tools.py
@@ -0,0 +1,15 @@
+""" Some grabbag tools """
+
+from pymbolic.primitives import Expression, Variable, Subscript
+
+
+def get_pymbolic_basename(expr):
+ assert isinstance(expr, Expression), "Type: {}, expr: {}".format(type(expr), expr)
+
+ if isinstance(expr, Variable):
+ return expr.name
+
+ if isinstance(expr, Subscript):
+ return get_pymbolic_basename(expr.aggregate)
+
+ raise NotImplementedError("Cannot determine basename of {}".format(expr))