diff --git a/python/dune/perftool/loopy/functions.py b/python/dune/perftool/loopy/functions.py deleted file mode 100644 index 703ffdc4f6cf0edb52926fb75a5509b73b77de1b..0000000000000000000000000000000000000000 --- a/python/dune/perftool/loopy/functions.py +++ /dev/null @@ -1,80 +0,0 @@ -from dune.perftool.generation import function_mangler -from loopy import CallMangleInfo -from loopy.symbolic import FunctionIdentifier -from loopy.types import NumpyType - -import numpy - - -class LFSChild(FunctionIdentifier): - def __init__(self, lfs): - self.lfs = lfs - - def __getinitargs__(self): - return (self.lfs,) - - @property - def name(self): - return '{}.child'.format(self.lfs) - - -@function_mangler -def lfs_child_mangler(target, func, dtypes): - if isinstance(func, LFSChild): - return CallMangleInfo(func.name, (NumpyType(str),), (NumpyType(numpy.int32),)) - - -class CoefficientAccess(FunctionIdentifier): - def __init__(self, container): - self.container = container - - def __getinitargs__(self): - return (self.container,) - - @property - def name(self): - return self.container - - -@function_mangler -def coefficient_mangler(target, func, dtypes): - if isinstance(func, CoefficientAccess): - return CallMangleInfo(func.name, (NumpyType(numpy.float64),), (NumpyType(str), NumpyType(numpy.int32))) - - -class PDELabAccumulationFunction(FunctionIdentifier): - def __init__(self, accumobj, rank): - self.accumobj = accumobj - self.rank = rank - - assert rank in (1, 2) - - def __getinitargs__(self): - return (self.accumobj, self.rank) - - @property - def name(self): - return '{}.accumulate'.format(self.accumobj) - - -@function_mangler -def accumulation_mangler(target, func, dtypes): - if isinstance(func, PDELabAccumulationFunction): - if func.rank == 1: - return CallMangleInfo(func.name, - (), - (NumpyType(str), - NumpyType(numpy.int32), - NumpyType(numpy.float64), - ) - ) - if func.rank == 2: - return CallMangleInfo(func.name, - (), - (NumpyType(str), - NumpyType(numpy.int32), - NumpyType(str), - NumpyType(numpy.int32), - NumpyType(numpy.float64), - ) - ) diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py index 02054f25d4d946b6096d690da3141785e5f1f8dd..d62f884dda2db26ce6d7f7c32c60c8136dc9b6d5 100644 --- a/python/dune/perftool/pdelab/argument.py +++ b/python/dune/perftool/pdelab/argument.py @@ -1,8 +1,13 @@ -""" Generator functions related to trial and test functions and the accumulation loop""" +""" Generator functions related to any input and output of accumulation kernels -from dune.perftool.options import get_option +Namely: +* Coefficient containers +* accumulation object (r, jac...) +""" +from dune.perftool.options import get_option from dune.perftool.generation import (domain, + function_mangler, iname, pymbolic_expr, globalarg, @@ -15,14 +20,76 @@ from dune.perftool.pdelab import (name_index, ) from dune.perftool.pdelab.basis import (evaluate_coefficient, evaluate_coefficient_gradient, - lfs_iname, name_basis, - name_lfs_bound, ) +from dune.perftool.pdelab.spaces import (lfs_iname, + name_lfs_bound, + ) from dune.perftool import Restriction + from pymbolic.primitives import Call, Subscript, Variable -import loopy +from loopy import CallMangleInfo +from loopy.symbolic import FunctionIdentifier +from loopy.types import NumpyType + +import numpy + + +class CoefficientAccess(FunctionIdentifier): + def __init__(self, container): + self.container = container + + def __getinitargs__(self): + return (self.container,) + + @property + def name(self): + return self.container + + +@function_mangler +def coefficient_mangler(target, func, dtypes): + if isinstance(func, CoefficientAccess): + return CallMangleInfo(func.name, (NumpyType(numpy.float64),), (NumpyType(str), NumpyType(numpy.int32))) + + +class PDELabAccumulationFunction(FunctionIdentifier): + def __init__(self, accumobj, rank): + self.accumobj = accumobj + self.rank = rank + + assert rank in (1, 2) + + def __getinitargs__(self): + return (self.accumobj, self.rank) + + @property + def name(self): + return '{}.accumulate'.format(self.accumobj) + + +@function_mangler +def accumulation_mangler(target, func, dtypes): + if isinstance(func, PDELabAccumulationFunction): + if func.rank == 1: + return CallMangleInfo(func.name, + (), + (NumpyType(str), + NumpyType(numpy.int32), + NumpyType(numpy.float64), + ) + ) + if func.rank == 2: + return CallMangleInfo(func.name, + (), + (NumpyType(str), + NumpyType(numpy.int32), + NumpyType(str), + NumpyType(numpy.int32), + NumpyType(numpy.float64), + ) + ) def name_trialfunction_gradient(element, restriction, component): @@ -67,22 +134,6 @@ def name_apply_function(element, restriction, component): return name -def name_testfunctionspace(restriction): - return restricted_name("lfsv", restriction) - - -def name_trialfunctionspace(restriction): - return restricted_name("lfsu", restriction) - - -def type_testfunctionspace(): - return "LFSV" - - -def type_trialfunctionspace(): - return "LFSU" - - def name_coefficientcontainer(restriction): name = restricted_name("x", restriction) return name @@ -97,14 +148,13 @@ def name_applycontainer(restriction): def pymbolic_coefficient(container, lfs, index): # TODO introduce a proper type for local function spaces! if isinstance(lfs, str): - valuearg(lfs, dtype=loopy.types.NumpyType("str")) + valuearg(lfs, dtype=NumpyType("str")) # If the LFS is not yet a pymbolic expression, make it one from pymbolic.primitives import Expression if not isinstance(lfs, Expression): lfs = Variable(lfs) - from dune.perftool.loopy.functions import CoefficientAccess return Call(CoefficientAccess(container), (lfs, Variable(index),)) @@ -112,21 +162,6 @@ def type_coefficientcontainer(): return "X" -def name_argumentspace(ma): - from ufl.classes import Argument, Coefficient - if isinstance(ma.argexpr, Argument): - if ma.argexpr.number() == 0: - return name_testfunctionspace(ma.restriction) - if ma.argexpr.number() == 1: - return name_trialfunctionspace(ma.restriction) - if isinstance(ma.argexpr, Coefficient): - # Index 0 is reserved for trialfunction, index 1 is reserved for jacobian apply function - assert ma.argexpr.count() < 2 - return name_trialfunctionspace(ma.restriction) - # 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): @@ -147,7 +182,6 @@ def type_residual(): def name_accumulation_variable(restrictions=(Restriction.NONE,)): - from dune.perftool.generation import get_global_context_value ft = get_global_context_value("form_type") if ft == 'residual' or ft == 'jacobian_apply': return name_residual(*restrictions) @@ -157,7 +191,6 @@ def name_accumulation_variable(restrictions=(Restriction.NONE,)): def type_accumulation_variable(): - from dune.perftool.generation import get_global_context_value ft = get_global_context_value("form_type") if ft == 'residual' or ft == 'jacobian_apply': return type_residual() diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py index 9d9e784feb81747f9ee7dbe9f06b5656eb49c303..b1ded62ea873994cfc521775d1a7e6da5693be07 100644 --- a/python/dune/perftool/pdelab/basis.py +++ b/python/dune/perftool/pdelab/basis.py @@ -2,14 +2,19 @@ from dune.perftool.generation import (cached, class_member, - domain, generator_factory, - iname, include_file, instruction, preamble, temporary_variable, ) +from dune.perftool.pdelab.spaces import (lfs_child, + lfs_iname, + name_leaf_lfs, + name_lfs, + name_lfs_bound, + type_gfs, + ) from dune.perftool.pdelab.quadrature import (name_quadrature_position_in_cell, quadrature_iname, ) @@ -22,6 +27,7 @@ from dune.perftool.pdelab.localoperator import (lop_template_ansatz_gfs, lop_template_test_gfs, ) from dune.perftool.pdelab.driver import FEM_name_mangling + from dune.perftool.pdelab import restricted_name from pymbolic.primitives import Product, Subscript, Variable from loopy import Reduction @@ -44,185 +50,6 @@ def name_localbasis_cache(element): return name -@preamble -def define_lfs_bound(lfs, bound): - return 'auto {} = {}.size();'.format(bound, lfs) - - -def name_lfs_bound(lfs): - # LFS might either be given by an UFL element or by a string describig its name - from ufl import FiniteElementBase - if isinstance(lfs, FiniteElementBase): - return name_lfs_bound(name_lfs(lfs)) - - bound = '{}_size'.format(lfs) - define_lfs_bound(lfs, bound) - return bound - - -@preamble -def using_indices(): - return "using namespace Dune::Indices;" - - -@preamble -def define_lfs(name, father, child): - using_indices() - return "auto {} = child({}, _{});".format(name, father, child) - - -def lfs_child(lfs, children, shape=None, symmetry=False): - from pymbolic.primitives import Call, Product, Sum - # Old pre-TensorElement implementation kept for comaptibility - if shape is None: - indices = (Variable(children[0]),) - else: - if symmetry and len(children) == 2: - # I do not want to think about tensors of rank > 2 right now - i, j = children - if i > j: - j, i = i, j - i = Variable(i) - j = Variable(j) - n = len(children) - - indices = (Sum((Product((n - 1, i)), Product((.5, i, 1 - i)), j)),) - else: - children = tuple(Variable(c) for c in children) - indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),) - - from dune.perftool.loopy.functions import LFSChild - return Call(LFSChild(lfs), indices) - - -@generator_factory(cache_key_generator=lambda e, r, **kw: (e, r)) -def name_leaf_lfs(leaf_element, restriction, val=None): - """ This function just caches leaf lfs names based on the - element. The resulting local function spaces are useful only - for size information. OTOH, they are available with just the - leaf element available (as seen in basis evaluation). - """ - # This generator function should be prepoluted by the lfs tree traversal, - # so val should always be None when we actually want the result - assert val - - return val - - -@generator_factory(cache_key_generator=lambda e, r, c, **kw: (e, r, c)) -def name_lfs(element, restriction, component, prefix=None): - # Omitting the prefix is only valid upon a second call, which will - # result in a cache hit. - assert prefix - - def _name_lfs(prefix, component): - name = prefix - if len(component) > 0: - name = name + '_' + '_'.join(str(i) for i in component) - return name - - name = _name_lfs(prefix, component) - if len(component) > 0: - father = _name_lfs(prefix, component[:-1]) - # If this localfunction space is the child of another one, trigger - # the extraction preamble. Necessary before going into recursion - # for having the correct (top-down) order of preambles - define_lfs(name, father, component[-1]) - - # Recurse into the given element to define all other local function spaces! - from ufl import MixedElement - from ufl.functionview import select_subelement - from ufl.classes import FixedIndex - subel = select_subelement(element, component) - if isinstance(subel, MixedElement): - for i in range(subel.num_sub_elements()): - name_lfs(element, restriction, component + (FixedIndex(i),), prefix=prefix) - - # Cache the name for the subelement - name_leaf_lfs(subel, restriction, val=name) - - # Now return the prefix! - return name - - -@generator_factory(cache_key_generator=lambda e, **kw: e) -def type_gfs(element, basetype=None, index_stack=None): - # Omitting basetype and index_stack is only valid upon a second call, - # which will result in a cache hit. - assert basetype - assert index_stack is not None - - # Additionally, element is expected to be a ufl finite element - from ufl import FiniteElementBase - assert isinstance(element, FiniteElementBase) - - # Recurse into the given element to define all other local function spaces! - from ufl import MixedElement - from ufl.classes import FixedIndex - if isinstance(element, MixedElement): - for i, subelem in enumerate(element.sub_elements()): - type_gfs(subelem, basetype=basetype, index_stack=index_stack + (FixedIndex(i),)) - - if len(index_stack) == 0: - return basetype - else: - include_file("dune/typetree/childextraction.hh", filetag="operatorfile") - return 'Dune::TypeTree::Child<{},{}>'.format(basetype, ','.join(str(i) for i in index_stack)) - - -def traverse_lfs_tree(arg): - from dune.perftool.ufl.modified_terminals import ModifiedArgumentDescriptor - assert isinstance(arg, ModifiedArgumentDescriptor) - - # First we need to determine the basename as given in the signature of - # this kernel method! - from dune.perftool.pdelab.argument import name_argumentspace - lfs_basename = name_argumentspace(arg) - from dune.perftool.pdelab.localoperator import lop_template_gfs - gfs_basename = lop_template_gfs(arg) - - # Now start recursively extracting local function spaces and fill the cache with - # all those values. That way we can later get a correct local function space with - # just the ufl finite element. - from ufl.classes import MultiIndex - name_lfs(arg.argexpr.ufl_element(), arg.restriction, MultiIndex(()), prefix=lfs_basename) - type_gfs(arg.argexpr.ufl_element(), basetype=gfs_basename, index_stack=()) - - -@generator_factory(item_tags=("iname",), cache_key_generator=lambda e, r, c: (e, c)) -def _lfs_iname(element, restriction, context): - lfs = name_leaf_lfs(element, restriction) - bound = name_lfs_bound(lfs) - - name = "{}_{}_index".format(lfs, context) - domain(name, bound) - - return name - - -def lfs_iname(element, restriction, count=None, context=''): - """ Get the iname to iterate over the local function space given by element - - Arguments: - ---------- - element: ufl.FiniteElementBase - The finite element this local function space belongs to - argcount: int - Use to realize double nesting in case of jacobians - context: str - Some generation methods will require you to duplicate an iname for - a given purpose, see the 'Loops and dependencies' of the loopy docs: - https://documen.tician.de/loopy/tutorial.html#loops-and-dependencies - """ - assert not ((context == '') and (count is None)) - if count is not None: - if context != '': - context = "{}_{}".format(count, context) - else: - context = str(count) - return _lfs_iname(element, restriction, context) - - @preamble def declare_cache_temporary(element, restriction, name, which): t_cache = type_localbasis_cache(element) diff --git a/python/dune/perftool/pdelab/signatures.py b/python/dune/perftool/pdelab/signatures.py index b3e09089a8dd3012412d26c2b025225f8210065a..3a40e636de6ba3224fdd1aaab2102419d3b4f1e6 100644 --- a/python/dune/perftool/pdelab/signatures.py +++ b/python/dune/perftool/pdelab/signatures.py @@ -9,11 +9,12 @@ from dune.perftool.pdelab.argument import (name_accumulation_variable, name_coefficientcontainer, type_coefficientcontainer, name_applycontainer, - name_testfunctionspace, - type_testfunctionspace, - name_trialfunctionspace, - type_trialfunctionspace, ) +from dune.perftool.pdelab.spaces import (name_testfunctionspace, + type_testfunctionspace, + name_trialfunctionspace, + type_trialfunctionspace, + ) def alpha_volume_signature(): diff --git a/python/dune/perftool/pdelab/spaces.py b/python/dune/perftool/pdelab/spaces.py new file mode 100644 index 0000000000000000000000000000000000000000..3812cc1fa49a52071748bc8788d1446fd6994166 --- /dev/null +++ b/python/dune/perftool/pdelab/spaces.py @@ -0,0 +1,243 @@ +""" Generator functions for PDELab local/grid function spaces etc. """ + +from dune.perftool.generation import (domain, + function_mangler, + generator_factory, + include_file, + preamble, + ) +from dune.perftool.pdelab import restricted_name + +from loopy import CallMangleInfo +from loopy.symbolic import FunctionIdentifier +from loopy.types import NumpyType + +from pymbolic.primitives import Variable + +import numpy + + +class LFSChild(FunctionIdentifier): + def __init__(self, lfs): + self.lfs = lfs + + def __getinitargs__(self): + return (self.lfs,) + + @property + def name(self): + return '{}.child'.format(self.lfs) + + +@function_mangler +def lfs_child_mangler(target, func, dtypes): + if isinstance(func, LFSChild): + return CallMangleInfo(func.name, (NumpyType(str),), (NumpyType(numpy.int32),)) + + +@preamble +def define_lfs_bound(lfs, bound): + return 'auto {} = {}.size();'.format(bound, lfs) + + +def name_lfs_bound(lfs): + # LFS might either be given by an UFL element or by a string describig its name + from ufl import FiniteElementBase + if isinstance(lfs, FiniteElementBase): + return name_lfs_bound(name_lfs(lfs)) + + bound = '{}_size'.format(lfs) + define_lfs_bound(lfs, bound) + return bound + + +@preamble +def using_indices(): + return "using namespace Dune::Indices;" + + +@preamble +def define_lfs(name, father, child): + using_indices() + return "auto {} = child({}, _{});".format(name, father, child) + + +def lfs_child(lfs, children, shape=None, symmetry=False): + from pymbolic.primitives import Call, Product, Sum + # Old pre-TensorElement implementation kept for comaptibility + if shape is None: + indices = (Variable(children[0]),) + else: + if symmetry and len(children) == 2: + # I do not want to think about tensors of rank > 2 right now + i, j = children + if i > j: + j, i = i, j + i = Variable(i) + j = Variable(j) + n = len(children) + + indices = (Sum((Product((n - 1, i)), Product((.5, i, 1 - i)), j)),) + else: + children = tuple(Variable(c) for c in children) + indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),) + + return Call(LFSChild(lfs), indices) + + +@generator_factory(cache_key_generator=lambda e, r, **kw: (e, r)) +def name_leaf_lfs(leaf_element, restriction, val=None): + """ This function just caches leaf lfs names based on the + element. The resulting local function spaces are useful only + for size information. OTOH, they are available with just the + leaf element available (as seen in basis evaluation). + """ + # This generator function should be prepoluted by the lfs tree traversal, + # so val should always be None when we actually want the result + assert val + + return val + + +@generator_factory(cache_key_generator=lambda e, r, c, **kw: (e, r, c)) +def name_lfs(element, restriction, component, prefix=None): + # Omitting the prefix is only valid upon a second call, which will + # result in a cache hit. + assert prefix + + def _name_lfs(prefix, component): + name = prefix + if len(component) > 0: + name = name + '_' + '_'.join(str(i) for i in component) + return name + + name = _name_lfs(prefix, component) + if len(component) > 0: + father = _name_lfs(prefix, component[:-1]) + # If this localfunction space is the child of another one, trigger + # the extraction preamble. Necessary before going into recursion + # for having the correct (top-down) order of preambles + define_lfs(name, father, component[-1]) + + # Recurse into the given element to define all other local function spaces! + from ufl import MixedElement + from ufl.functionview import select_subelement + from ufl.classes import FixedIndex + subel = select_subelement(element, component) + if isinstance(subel, MixedElement): + for i in range(subel.num_sub_elements()): + name_lfs(element, restriction, component + (FixedIndex(i),), prefix=prefix) + + # Cache the name for the subelement + name_leaf_lfs(subel, restriction, val=name) + + # Now return the prefix! + return name + + +@generator_factory(cache_key_generator=lambda e, **kw: e) +def type_gfs(element, basetype=None, index_stack=None): + # Omitting basetype and index_stack is only valid upon a second call, + # which will result in a cache hit. + assert basetype + assert index_stack is not None + + # Additionally, element is expected to be a ufl finite element + from ufl import FiniteElementBase + assert isinstance(element, FiniteElementBase) + + # Recurse into the given element to define all other local function spaces! + from ufl import MixedElement + from ufl.classes import FixedIndex + if isinstance(element, MixedElement): + for i, subelem in enumerate(element.sub_elements()): + type_gfs(subelem, basetype=basetype, index_stack=index_stack + (FixedIndex(i),)) + + if len(index_stack) == 0: + return basetype + else: + include_file("dune/typetree/childextraction.hh", filetag="operatorfile") + return 'Dune::TypeTree::Child<{},{}>'.format(basetype, ','.join(str(i) for i in index_stack)) + + +def traverse_lfs_tree(arg): + from dune.perftool.ufl.modified_terminals import ModifiedArgumentDescriptor + assert isinstance(arg, ModifiedArgumentDescriptor) + + # First we need to determine the basename as given in the signature of + # this kernel method! + lfs_basename = name_argumentspace(arg) + from dune.perftool.pdelab.localoperator import lop_template_gfs + gfs_basename = lop_template_gfs(arg) + + # Now start recursively extracting local function spaces and fill the cache with + # all those values. That way we can later get a correct local function space with + # just the ufl finite element. + from ufl.classes import MultiIndex + name_lfs(arg.argexpr.ufl_element(), arg.restriction, MultiIndex(()), prefix=lfs_basename) + type_gfs(arg.argexpr.ufl_element(), basetype=gfs_basename, index_stack=()) + + +@generator_factory(item_tags=("iname",), cache_key_generator=lambda e, r, c: (e, c)) +def _lfs_iname(element, restriction, context): + lfs = name_leaf_lfs(element, restriction) + bound = name_lfs_bound(lfs) + + name = "{}_{}_index".format(lfs, context) + domain(name, bound) + + return name + + +def lfs_iname(element, restriction, count=None, context=''): + """ Get the iname to iterate over the local function space given by element + + Arguments: + ---------- + element: ufl.FiniteElementBase + The finite element this local function space belongs to + argcount: int + Use to realize double nesting in case of jacobians + context: str + Some generation methods will require you to duplicate an iname for + a given purpose, see the 'Loops and dependencies' of the loopy docs: + https://documen.tician.de/loopy/tutorial.html#loops-and-dependencies + """ + assert not ((context == '') and (count is None)) + if count is not None: + if context != '': + context = "{}_{}".format(count, context) + else: + context = str(count) + return _lfs_iname(element, restriction, context) + + +def name_testfunctionspace(restriction): + return restricted_name("lfsv", restriction) + + +def name_trialfunctionspace(restriction): + return restricted_name("lfsu", restriction) + + +def name_argumentspace(ma): + from ufl.classes import Argument, Coefficient + if isinstance(ma.argexpr, Argument): + if ma.argexpr.number() == 0: + return name_testfunctionspace(ma.restriction) + if ma.argexpr.number() == 1: + return name_trialfunctionspace(ma.restriction) + if isinstance(ma.argexpr, Coefficient): + # Index 0 is reserved for trialfunction, index 1 is reserved for jacobian apply function + assert ma.argexpr.count() < 2 + return name_trialfunctionspace(ma.restriction) + # We should never encounter an argument other than 0 or 1 + assert False + + +def type_testfunctionspace(): + return "LFSV" + + +def type_trialfunctionspace(): + return "LFSU" diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py index cfc03a1f0714b44582f7ca041812f446480065fe..26c381b6bc4029a84e80784f2557ff137d981f2c 100644 --- a/python/dune/perftool/ufl/visitor.py +++ b/python/dune/perftool/ufl/visitor.py @@ -16,16 +16,17 @@ from dune.perftool.generation import (domain, valuearg, ) -from dune.perftool.pdelab.basis import (lfs_iname, - name_leaf_lfs, - name_lfs, - name_lfs_bound, - traverse_lfs_tree, - ) +from dune.perftool.pdelab.spaces import (lfs_iname, + name_leaf_lfs, + name_lfs, + name_lfs_bound, + traverse_lfs_tree, + ) from dune.perftool.pdelab.quadrature import quadrature_iname from pymbolic.primitives import Subscript, Variable from ufl.algorithms import MultiFunction + class UFL2LoopyVisitor(ModifiedTerminalTracker, GeometryMapper): def __init__(self, measure, subdomain_id, dimension_indices): # Some variables describing the integral measure of this integral @@ -59,6 +60,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, GeometryMapper): if self.measure == 'exterior_facet': ma.restriction = Restriction.NEGATIVE + from dune.perftool.pdelab.spaces import traverse_lfs_tree traverse_lfs_tree(ma) # Determine the rank of the term @@ -204,6 +206,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, GeometryMapper): raise ValueError("Gradients should have been transformed to reference gradients!!!") # Have the issued instruction depend on the iname for this localfunction space + from dune.perftool.pdelab.spaces import lfs_iname iname = lfs_iname(leaf_element, restriction, o.number()) self.inames.append(iname)