localoperator.py

from __future__ import absolute_import

from dune.perftool.options import get_option
from dune.perftool.generation import (base_class,
                                      class_basename,
                                      class_member,
                                      constructor_parameter,
                                      global_context,
                                      include_file,
                                      initializer_list,
                                      symbol,
                                      template_parameter,
                                      )
from dune.perftool.cgen.clazz import (AccessModifier,
                                      BaseClass,
                                      ClassMember,
                                      )
from dune.perftool import Restriction


@template_parameter("operator")
def lop_template_ansatz_gfs():
    return "GFSU"


@template_parameter("operator")
def lop_template_test_gfs():
    return "GFSV"


def lop_template_gfs(ma):
    from ufl.classes import Argument, Coefficient
    if isinstance(ma.argexpr, Argument):
        if ma.argexpr.number() == 0:
            return lop_template_test_gfs()
        if ma.argexpr.number() == 1:
            return lop_template_ansatz_gfs()
    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 lop_template_ansatz_gfs()
    assert False


@symbol
def name_initree_constructor_param():
    return "iniParams"


@class_member("operator")
def define_initree(name):
    param_name = name_initree_constructor_param()
    include_file('dune/common/parametertree.hh', filetag="operatorfile")
    constructor_parameter("const Dune::ParameterTree&", param_name, classtag="operator")
    initializer_list(name, [param_name], classtag="operator")

    return "const Dune::ParameterTree& {};".format(name)


def ufl_measure_to_pdelab_measure(which):
    return {'cell': 'Volume',
            'exterior_facet': 'Boundary',
            'interior_facet': 'Skeleton',
            }.get(which)


@class_member(classtag="operator", access=AccessModifier.PUBLIC)
def _enum_pattern(which):
    return "enum {{ doPattern{} = true }};".format(which)


def enum_pattern():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    return _enum_pattern(ufl_measure_to_pdelab_measure(integral_type))


def _pattern_baseclass(measure):
    return base_class('Dune::PDELab::Full{}Pattern'.format(measure), classtag="operator")


def pattern_baseclass():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    return _pattern_baseclass(ufl_measure_to_pdelab_measure(integral_type))


@class_member(classtag="operator", access=AccessModifier.PUBLIC)
def _enum_alpha(which):
    return "enum {{ doAlpha{} = true }};".format(which)


def enum_alpha():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    return _enum_alpha(ufl_measure_to_pdelab_measure(integral_type))


@symbol
def name_initree_member():
    define_initree("_iniParams")
    return "_iniParams"


@class_basename("operator")
def localoperator_basename():
    return "LocalOperator"


def class_type_from_cache(classtag):
    from dune.perftool.generation import retrieve_cache_items

    # get the basename
    basename = [i for i in retrieve_cache_items(condition="{} and basename".format(classtag))]
    assert len(basename) == 1
    basename = basename[0]

    # get the template parameters
    tparams = [i for i in retrieve_cache_items(condition="{} and template_param".format(classtag))]
    tparam_str = ''
    if len(tparams) > 0:
        tparam_str = '<{}>'.format(', '.join(t for t in tparams))

    return basename, basename + tparam_str


def assembly_routine_signature():
    from dune.perftool.generation import get_global_context_value
    integral_type = get_global_context_value("integral_type")
    form_type = get_global_context_value("form_type")

    if form_type == 'residual':
        if integral_type == 'cell':
            from dune.perftool.pdelab.signatures import alpha_volume_signature
            return alpha_volume_signature()
        if integral_type == 'exterior_facet':
            from dune.perftool.pdelab.signatures import alpha_boundary_signature
            return alpha_boundary_signature()
        if integral_type == 'interior_facet':
            from dune.perftool.pdelab.signatures import alpha_skeleton_signature
            return alpha_skeleton_signature()

    if form_type == 'jacobian':
        if integral_type == 'cell':
            from dune.perftool.pdelab.signatures import jacobian_volume_signature
            return jacobian_volume_signature()
        if integral_type == 'exterior_facet':
            from dune.perftool.pdelab.signatures import jacobian_boundary_signature
            return jacobian_boundary_signature()
        if integral_type == 'interior_facet':
            from dune.perftool.pdelab.signatures import jacobian_skeleton_signature
            return jacobian_skeleton_signature()

    if form_type == 'jacobian_apply':
        if integral_type == 'cell':
            from dune.perftool.pdelab.signatures import jacobian_apply_volume_signature
            return jacobian_apply_volume_signature()
        if integral_type == 'exterior_facet':
            from dune.perftool.pdelab.signatures import jacobian_apply_boundary_signature
            return jacobian_apply_boundary_signature()
        if integral_type == 'interior_facet':
            from dune.perftool.pdelab.signatures import jacobian_apply_skeleton_signature
            return jacobian_apply_skeleton_signature()

    assert False


def generate_kernel(integrals):
    subst_rules = []
    for integral in integrals:
        integrand = integral.integrand()
        measure = integral.integral_type()
        subdomain_id = integral.subdomain_id()
        subdomain_data = integral.subdomain_data()

        from dune.perftool.ufl.dimensionindex import dimension_index_mapping
        dimension_indices = dimension_index_mapping(integrand)

        # Now split the given integrand into accumulation expressions
        from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms
        accterms = split_into_accumulation_terms(integrand)

        # Get a transformer instance for this kernel
        from dune.perftool.loopy.transformer import UFL2LoopyVisitor
        visitor = UFL2LoopyVisitor(measure, subdomain_id, dimension_indices)

        # Iterate over the terms and generate a kernel
        for term in accterms:
            visitor(term)
            subst_rules.extend(visitor.substitution_rules)

    # Extract the information, which is needed to create a loopy kernel.
    # First extracting it, might be useful to alter it before kernel generation.
    from dune.perftool.generation import retrieve_cache_items
    from dune.perftool.loopy.target import DuneTarget
    domains = [i for i in retrieve_cache_items("domain")]
    instructions = [i for i in retrieve_cache_items("instruction")]
    temporaries = {i.name: i for i in retrieve_cache_items("temporary")}
    arguments = [i for i in retrieve_cache_items("argument")]

    # Get the function manglers
    from dune.perftool.loopy.functions import accumulation_mangler, coefficient_mangler, lfs_child_mangler

    # Create the kernel
    from loopy import make_kernel, preprocess_kernel
    kernel = make_kernel(domains,
                         instructions + subst_rules,
                         arguments,
                         temporary_variables=temporaries,
                         function_manglers=[accumulation_mangler, coefficient_mangler, lfs_child_mangler],
                         target=DuneTarget()
                         )

    from loopy import make_reduction_inames_unique
    kernel = make_reduction_inames_unique(kernel)

    kernel = preprocess_kernel(kernel)

    # see whether we need iname duplication to make this thing schedulable!
    # TODO: Use a clever strategy here, instead of random transformation until the problem is resolved
    from loopy import needs_iname_duplication, get_iname_duplication_options, duplicate_inames
    while needs_iname_duplication(kernel):
        # If there is a duplication that solves the problem with just one duplication, we pick that one
        iname, within = (None, None)
        for i, w in get_iname_duplication_options(kernel):
            if not needs_iname_duplication(duplicate_inames(kernel, i, w)):
                iname, within = (i, w)

        # Otherwise pick a random one.
        if iname is None:
            iname, within = next(get_iname_duplication_options(kernel))

        # Do the transformation
        print "Applying iname duplication to measure {}: iname {}; within {}".format(measure, iname, within)
        kernel = duplicate_inames(kernel, iname, within)

    # Loopy might have introduced some temporary variables during preprocessing. As I want to have my own
    # temporary declaration code right now, I call the declaration preamble manually.
    for added_tv in set(kernel.temporary_variables.keys()) - set(temporaries.keys()):
        from dune.perftool.generation.loopy import default_declaration
        default_declaration(added_tv)

    # Now add the preambles to the kernel
    preambles = [(i, p) for i, p in enumerate(retrieve_cache_items("preamble"))]
    kernel = kernel.copy(preambles=preambles)

    # All items with the kernel tags can be destroyed once a kernel has been generated
    from dune.perftool.generation import delete_cache_items
    delete_cache_items("(not file) and (not clazz)")

    return kernel


class AssemblyMethod(ClassMember):
    def __init__(self, signature, kernel):
        from loopy import generate_body
        from cgen import LiteralLines, Block
        content = signature
        content.append('{')
        if kernel is not None:
            for i, p in kernel.preambles:
                content.append('  ' + p)
            content.extend(l for l in generate_body(kernel).split('\n')[1:-1])
        content.append('}')
        ClassMember.__init__(self, content)


def cgen_class_from_cache(tag, members=[]):
    from dune.perftool.generation import retrieve_cache_items

    # Generate the name by concatenating basename and template parameters
    basename, fullname = class_type_from_cache(tag)

    base_classes = [bc for bc in retrieve_cache_items('{} and baseclass'.format(tag))]
    constructor_params = [bc for bc in retrieve_cache_items('{} and constructor_param'.format(tag))]
    il = [i for i in retrieve_cache_items('{} and initializer'.format(tag))]
    pm = [m for m in retrieve_cache_items('{} and member'.format(tag))]
    tparams = [i for i in retrieve_cache_items('{} and template_param'.format(tag))]

    from dune.perftool.cgen.clazz import Constructor
    constructor = Constructor(arg_decls=constructor_params, clsname=basename, initializer_list=il)

    from dune.perftool.cgen import Class
    return Class(basename, base_classes=base_classes, members=members + pm, constructors=[constructor], tparam_decls=tparams)


def generate_localoperator_kernels(formdata, data):
    # Extract the relevant attributes of the form data
    form = formdata.preprocessed_form

    # Reset the generation cache
    from dune.perftool.generation import delete_cache_items
    delete_cache_items()

    # Manage includes and base classes that we always need
    include_file('dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh', filetag="operatorfile")
    include_file('dune/pdelab/localoperator/idefault.hh', filetag="operatorfile")
    include_file('dune/pdelab/localoperator/flags.hh', filetag="operatorfile")
    include_file('dune/pdelab/localoperator/pattern.hh', filetag="operatorfile")

    # Trigger this one once early on to avoid wrong stuff happening
    localoperator_basename()
    lop_template_ansatz_gfs()
    lop_template_test_gfs()
    from dune.perftool.pdelab.parameter import parameterclass_basename
    parameterclass_basename()

    # Make sure there is always the same constructor arguments (even if parameter class is empty)
    from dune.perftool.pdelab.localoperator import name_initree_member
    name_initree_member()
    from dune.perftool.pdelab.parameter import name_paramclass
    name_paramclass()

    base_class('Dune::PDELab::LocalOperatorDefaultFlags', classtag="operator")

    # Have a data structure collect the generated kernels
    operator_kernels = {}

    with global_context(form_type='residual'):
        # Generate the necessary residual methods
        for measure in set(i.integral_type() for i in form.integrals()):
            with global_context(integral_type=measure):
                enum_pattern()
                pattern_baseclass()
                enum_alpha()
                kernel = generate_kernel(form.integrals_by_type(measure))

                # Maybe add numerical differentiation
                if get_option("numerical_jacobian"):
                    include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")

                    _, loptype = class_type_from_cache("operator")
                    which = ufl_measure_to_pdelab_measure(measure)

                    # Numerical jacobian methods
                    base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")

                    # Add the initializer list for that base class
                    ini = name_initree_member()
                    ini_constructor = name_initree_constructor_param()

                    # Numerical jacobian methods
                    initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
                                     ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
                                     classtag="operator",
                                     )

                    if get_option("matrix_free"):
                        # Numeical jacobian apply base class
                        base_class("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype), classtag="operator")

                        # Numerical jacobian apply initializer list
                        initializer_list("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype),
                                         ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
                                         classtag="operator",
                                         )

                operator_kernels[(measure, 'residual')] = kernel

    # 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]))

        with global_context(form_type="jacobian"):
            for measure in set(i.integral_type() for i in jacform.integrals()):
                with global_context(integral_type=measure):
                    kernel = generate_kernel(jacform.integrals_by_type(measure))
                operator_kernels[(measure, 'jacobian')] = kernel

            # Generate dummy functions for those kernels, that vanished in the differentiation process
            # We *could* solve this problem by using lambda_* terms but we do not really want that, so
            # we use empty jacobian assembly methods instead
            alpha_measures = set(i.integral_type() for i in form.integrals())
            jacobian_measures = set(i.integral_type() for i in jacform.integrals())
            for it in alpha_measures - jacobian_measures:
                operator_kernels[(it, 'jacobian')] = None

        # Jacobian apply methods for matrix-free computations
        if get_option("matrix_free"):
            # The apply vector has reserved index 1 so we directly use Coefficient class from ufl
            from ufl import Coefficient
            apply_coefficient = Coefficient(form.coefficients()[0].ufl_element(), 1)

            # Create application of jacobian on vector
            from ufl import action
            jac_apply_form = action(jacform, apply_coefficient)

            # Create kernel for jacobian application
            with global_context(form_type="jacobian_apply"):
                for measure in set(i.integral_type() for i in jac_apply_form.integrals()):
                    with global_context(integral_type=measure):
                        kernel = generate_kernel(jac_apply_form.integrals_by_type(measure))
                    operator_kernels[(measure, 'jacobian_apply')] = kernel

                    # Generate dummy functions for those kernels, that vanished in the differentiation process
                    # We *could* solve this problem by using lambda_* terms but we do not really want that, so
                    # we use empty jacobian assembly methods instead
                    alpha_measures = set(i.integral_type() for i in form.integrals())
                    jacobian_apply_measures = set(i.integral_type() for i in jac_apply_form.integrals())
                    for it in alpha_measures - jacobian_apply_measures:
                        operator_kernels[(it, 'jacobian_apply')] = None

    # Return the set of generated kernels
    return operator_kernels


def generate_localoperator_file(kernels):
    operator_methods = []

    # Make generables from the given kernels
    for method, kernel in kernels.items():
        it, ft = method
        with global_context(integral_type=it, form_type=ft):
            signature = assembly_routine_signature()
            operator_methods.append(AssemblyMethod(signature, kernel))

    # Write the file!
    from dune.perftool.file import generate_file

    param = cgen_class_from_cache("parameterclass")
    # TODO take the name of this thing from the UFL file
    lop = cgen_class_from_cache("operator", members=operator_methods)
    generate_file(get_option("operator_file"), "operatorfile", [param, lop])