localoperator.py

from pytools import memoize
from dune.perftool.options import get_option
from dune.perftool.generation import generator_factory
from dune.perftool.pdelab import dune_symbol

# Define the generators used in-here
operator_include = generator_factory(item_tags=("pdelab", "include", "operator"), on_store=lambda i: "#include<{}>".format(i), no_deco=True)
base_class = generator_factory(item_tags=("pdelab", "baseclass", "operator"), counted=True, no_deco=True)
initializer_list = generator_factory(item_tags=("pdelab", "initializer", "operator"), counted=True, no_deco=True)
# TODO definition
#private_member = generator_factory(item_tags=("pdelab", "member", "privatemember"))

@generator_factory(item_tags=("pdelab", "constructor"), counted=True)
def constructor_parameter(_type, name):
    return "{} {}".format(_type, name)

@dune_symbol
def name_initree_constructor():
    operator_include('dune/common/parametertree.hh')
    constructor_parameter("const Dune::ParameterTree&", "iniParams")
    return "iniParams"

@memoize
def measure_specific_details(measure):
    # The return dictionary that this memoized method will grant direct access to.
    ret = {}

    def numerical_jacobian(which):
        if get_option("numerical_jacobian"):
            # Add a base class
            from dune.perftool.pdelab.driver import type_localoperator
            loptype = type_localoperator()
            base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype))

            # Add the initializer list for that base class
            ini = name_initree_constructor()
            initializer_list("Dune::PDELab::NumericalJacobian{}<{}>({}.get(\"numerical_epsilon.{}\", 1e-9))".format(which, loptype, ini, which.lower()))


    if measure == "cell":
        base_class('Dune::PDELab::FullVolumePattern')
        numerical_jacobian("Volume")

        ret["residual_signature"] = ['template<typename EG, typename LFSV0, typename X, typename LFSV1, typename R>',
                                     'void alpha_volume(const EG& eg, const LFSV0& lfsv0, const X& x, const LFSV1& lfsv1, R& r) const']
        ret["jacobian_signature"] = ['template<typename EG, typename LFSV0, typename X, typename LFSV1, typename J>',
                                     'void jacobian_volume(const EG& eg, const LFSV0& lfsv0, const X& x, const LFSV1& lfsv1, J& jac) const']

    if measure == "exterior_facet":
        base_class('Dune::PDELab::FullBoundaryPattern')
        numerical_jacobian("Boundary")

        ret["residual_signature"] = ['template<typename IG, typename LFSV0, typename X, typename LFSV1, typename R>',
                                     'void alpha_boundary(const IG& ig, const LFSV0& lfsv0, const X& x, const LFSV1& lfsv1, R& r) const']
        ret["jacobian_signature"] = ['template<typename IG, typename LFSV0, typename X, typename LFSV1, typename J>',
                                     'void jacobian_boundary(const IG& ig, const LFSV0& lfsv0, const X& x, const LFSV1& lfsv1, J& jac) const']

    if measure == "interior_facet":
        base_class('Dune::PDELab::FullSkeletonPattern')
        numerical_jacobian("Skeleton")

        ret["residual_signature"] = ['template<typename IG, typename LFSV0_S, typename X, typename LFSV1_S, typename LFSV0_N, typename R, typename LFSV1_N>',
                                     'void alpha_skeleton(const IG& ig, const LFSV0_S& lfsv0_s, const X& x_s, const LFSV1_S& lfsv1_s, const LFSV0_N& lfsv0_n, const X& x_n, const LFSV1_N& lfsv1_n, R& r_s, R& r_n) const']
        ret["jacobian_signature"] = ['template<typename IG, typename LFSV0_S, typename X, typename LFSV1_S, typename LFSV0_N, typename LFSV1_N, typename Jac>',
                                     'void jacobian_skeleton(const IG& ig, const LFSV0_S& lfsv0_s, const X& x_s, const LFSV1_S& lfsv1_s, const LFSV0_N& lfsv0_n, const X& x_n, const LFSV1_N& lfsv1_n, Jac& jac_ss, Jac& jac_sn, Jac& jac_ns, Jac& jac_nn) const']

    return ret


def generate_term(integrand=None, measure=None):
    assert integrand and measure

    # Delete all non-include parts of the cache.
    # TODO: add things such as base classes as cache items.
    from dune.perftool.generation import delete_cache
    delete_cache()

    # Get the measure specifics
    specifics = measure_specific_details(measure)

    # Transform the expression. All relevant results are then in the cache
    from dune.perftool.transformer import transform_expression
    transform_expression(integrand)

    # 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.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")}
    preambles = [i for i in retrieve_cache_items("preamble")]

# TODO it might be necessary to list the arguments and their shape manually to avoid automatic shape detection!
    import loopy
    import numpy
    from loopy import make_kernel, preprocess_kernel
    from dune.perftool.target import DuneTarget
    k = make_kernel(domains, instructions,
                [loopy.GlobalArg("grad_a0", numpy.float64, shape=("argi_n", "dim")),
                 loopy.GlobalArg("grad_a1", numpy.float64, shape=("argj_n", "dim")),
                 loopy.GlobalArg("grad_c0", numpy.float64, shape=("dim")),
                 loopy.ValueArg("dim", numpy.int32),
                 loopy.ValueArg("q_n", numpy.int32),
                 loopy.ValueArg("argi_n", numpy.int32),
                 loopy.ValueArg("argj_n", numpy.int32)],
                    temporary_variables=temporaries, target=DuneTarget())
    k = preprocess_kernel(k)
    kernel = k

#     # Create the kernel
#     from loopy import make_kernel, preprocess_kernel, add_dtypes
#     kernel = make_kernel(domains, instructions, temporary_variables=temporaries, preambles=preambles, target=DuneTarget())
#     import numpy
#     kernel = add_dtypes(kernel, dict(grad_a0=numpy.float64, grad_c0=numpy.float64))
#     kernel = preprocess_kernel(kernel)

    # Return the actual code (might instead return kernels...)
    from loopy import generate_code
    return str(generate_code(kernel)[0])


def generate_localoperator(ufldata, operatorfile):
    form = ufldata.preprocessed_form

    operator_methods = []

    # For the moment, I do assume that there is but one integral of each type. This might differ
    # if you use different quadrature orders for different terms.
    assert len(form.integrals()) == len(set(i.integral_type() for i in form.integrals()))

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

    # Generate the necessary residual methods
    for integral in form.integrals():
        body = generate_term(integrand=integral.integrand(), measure=integral.integral_type())
        signature = measure_specific_details(integral.integral_type())["residual_signature"]
        operator_methods.append((signature, body))

    # Generate the necessary jacobian methods
    from dune.perftool.options import get_option
    if get_option("numerical_jacobian"):
        operator_include("dune/pdelab/localoperator/defaultimp.hh")
    else:
        from ufl import derivative
        from ufl.algorithms import expand_derivatives
        jacform = expand_derivatives(derivative(form, form.coefficients()[0]))

        for integral in jacform.integrals():
            body = generate_term(integrand=integral.integrand(), measure=integral.integral_type())
            signature = measure_specific_details(integral.integral_type())["jacobian_signature"]
            operator_methods.append((signature, body))


    # TODO: JacobianApply for matrix-free computations.

    # Manage includes and base classes that we always need
    operator_include('dune/pdelab/gridfunctionspace/gridfunctionspaceutilities.hh')
    operator_include('dune/pdelab/localoperator/idefault.hh')
    operator_include('dune/pdelab/localoperator/flags.hh')
    operator_include('dune/pdelab/localoperator/pattern.hh')
    operator_include('dune/geometry/quadraturerules.hh')

    base_class('Dune::PDELab::LocalOperatorDefaultFlags')