diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py index 28019f5845a1473f6bd40ff5594a0c0b9040ebc3..3301357bffddb1a4d1815e87e240f01d00cce14b 100644 --- a/python/dune/perftool/pdelab/localoperator.py +++ b/python/dune/perftool/pdelab/localoperator.py @@ -5,6 +5,7 @@ from dune.perftool.generation import (base_class, class_basename, class_member, constructor_parameter, + global_context, include_file, initializer_list, symbol, @@ -43,16 +44,45 @@ def define_initree(name): 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): +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): +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") @@ -81,58 +111,26 @@ def class_type_from_cache(classtag): return basename, basename + tparam_str -@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 - _, loptype = class_type_from_cache("operator") - 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() - initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), - ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())], - classtag="operator") - - if measure == "cell": - base_class('Dune::PDELab::FullVolumePattern', classtag="operator") - numerical_jacobian("Volume") - enum_pattern("Volume") - enum_alpha("Volume") - - ret["residual_signature"] = ['template<typename EG, typename LFSU, typename X, typename LFSV, typename R>', - 'void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const'] - ret["jacobian_signature"] = ['template<typename EG, typename LFSU, typename X, typename LFSV, typename J>', - 'void jacobian_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, J& jac) const'] - - if measure == "exterior_facet": - base_class('Dune::PDELab::FullBoundaryPattern', classtag="operator") - numerical_jacobian("Boundary") - enum_pattern("Boundary") - enum_alpha("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', classtag="operator") - numerical_jacobian("Skeleton") - enum_pattern("Skeleton") - enum_alpha("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 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") + + aj = "alpha" if form_type == 'residual' else "jacobian" + + if integral_type == 'cell': + return ['template<typename EG, typename LFSU, typename X, typename LFSV, typename R>', + 'void {}_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const'.format(aj)] + + if integral_type == 'exterior_facet': + return ['template<typename IG, typename LFSV0, typename X, typename LFSV1, typename R>', + 'void {}_boundary(const IG& ig, const LFSV0& lfsv0, const X& x, const LFSV1& lfsv1, R& r) const'.format(aj)] + + if integral_type == 'interior_facet': + return ['template<typename IG, typename LFSV0_S, typename X, typename LFSV1_S, typename LFSV0_N, typename R, typename LFSV1_N>', + 'void {}_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'.format(aj)] + + assert False def generate_kernel(integral): @@ -141,9 +139,6 @@ def generate_kernel(integral): subdomain_id = integral.subdomain_id() subdomain_data = integral.subdomain_data() - # Get the measure specifics - specifics = measure_specific_details(measure) - # 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) @@ -245,23 +240,44 @@ def generate_localoperator_kernels(formdata, namedata): from dune.perftool.generation import global_context with global_context(formdata=formdata, namedata=namedata): - # Generate the necessary residual methods - for integral in form.integrals(): - kernel = generate_kernel(integral) - operator_kernels[(integral.integral_type(), 'residual')] = kernel + with global_context(form_type='residual'): + # Generate the necessary residual methods + for integral in form.integrals(): + with global_context(integral_type=integral.integral_type()): + enum_pattern() + pattern_baseclass() + enum_alpha() + kernel = generate_kernel(integral) + + # 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(integral.integral_type()) + 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() + initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), + ["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())], + classtag="operator", + ) + + operator_kernels[(integral.integral_type(), 'residual')] = kernel # Generate the necessary jacobian methods - from dune.perftool.options import get_option - if get_option("numerical_jacobian"): - include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile") - else: + if not get_option("numerical_jacobian"): from ufl import derivative from ufl.algorithms import expand_derivatives jacform = expand_derivatives(derivative(form, form.coefficients()[0])) - for integral in jacform.integrals(): - kernel = generate_kernel(integral) - operator_kernels[(integral.integral_type(), 'jacobian')] = kernel + with global_context(form_type="jacobian"): + for integral in jacform.integrals(): + with global_context(integral_type=integral.integral_type()): + kernel = generate_kernel(integral) + operator_kernels[(integral.integral_type(), 'jacobian')] = kernel # TODO: JacobianApply for matrix-free computations. @@ -274,8 +290,10 @@ def generate_localoperator_file(kernels): # Make generables from the given kernels for method, kernel in kernels.items(): - signature = measure_specific_details(method[0])["{}_signature".format(method[1])] - operator_methods.append(AssemblyMethod(signature, kernel)) + 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