diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py index f948414ab1728b8266a35f0d2ea48853f5e487c0..5a044e5532f887d82ad646959c98a59f24a6c916 100644 --- a/python/dune/perftool/pdelab/argument.py +++ b/python/dune/perftool/pdelab/argument.py @@ -176,14 +176,14 @@ def name_accumulation_variable(restrictions=None): if measure == "cell": restrictions = (Restriction.NONE,) else: - restrictions = (Restriction.OUTSIDE,) + restrictions = (Restriction.NEGATIVE,) return name_residual(*restrictions) if ft == 'jacobian': if restrictions is None: if measure == "cell": restrictions = (Restriction.NONE, Restriction.NONE) else: - restrictions = (Restriction.OUTSIDE, Restriction.OUTSIDE) + restrictions = (Restriction.NEGATIVE, Restriction.NEGATIVE) return name_jacobian(*restrictions) assert False diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py index d39bd9cb578b62e17f64eb6920bcde855cc32980..3f8392af0939fcda445137c487d7fbf4ef17e10b 100644 --- a/python/dune/perftool/pdelab/geometry.py +++ b/python/dune/perftool/pdelab/geometry.py @@ -232,6 +232,10 @@ def name_dimension(): return formdata.geometric_dimension +def world_dimension(): + return name_dimension() + + def name_intersection_dimension(): formdata = get_global_context_value('formdata') return formdata.geometric_dimension - 1 diff --git a/python/dune/perftool/pdelab/signatures.py b/python/dune/perftool/pdelab/signatures.py index d903f929436fd085ea3ed7012f09b77c32fd7568..d2c48409d746278ea35d9abb22d99a9bad6675e8 100644 --- a/python/dune/perftool/pdelab/signatures.py +++ b/python/dune/perftool/pdelab/signatures.py @@ -56,12 +56,12 @@ def assembly_routine_signature(): formdata = get_global_context_value("formdata") templates, args = {('residual', 'cell'): (alpha_volume_templates, alpha_volume_args), - ('residual', 'exterior_facet'): (alpha_boundary_templates, alpha_boundary_args), - ('residual', 'interior_facet'): (alpha_skeleton_templates, alpha_skeleton_args), - ('jacobian', 'cell'): (jacobian_volume_templates, jacobian_volume_args), - ('jacobian', 'exterior_facet'): (jacobian_boundary_templates, jacobian_boundary_args), - ('jacobian', 'interior_facet'): (jacobian_skeleton_templates, jacobian_skeleton_args), - }.get((form_type, integral_type), (None, None)) + ('residual', 'exterior_facet'): (alpha_boundary_templates, alpha_boundary_args), + ('residual', 'interior_facet'): (alpha_skeleton_templates, alpha_skeleton_args), + ('jacobian', 'cell'): (jacobian_volume_templates, jacobian_volume_args), + ('jacobian', 'exterior_facet'): (jacobian_boundary_templates, jacobian_boundary_args), + ('jacobian', 'interior_facet'): (jacobian_skeleton_templates, jacobian_skeleton_args), + }.get((form_type, integral_type), (None, None)) if templates is None: # Check if form is linear @@ -69,16 +69,45 @@ def assembly_routine_signature(): linear = is_linear(formdata.original_form) templates, args = {('jacobian_apply', 'cell', True): (jacobian_apply_volume_templates, jacobian_apply_volume_args), - ('jacobian_apply', 'exterior_facet', True): (jacobian_apply_boundary_templates, jacobian_apply_boundary_args), - ('jacobian_apply', 'interior_facet', True): (jacobian_apply_skeleton_templates, jacobian_apply_skeleton_args), - ('jacobian_apply', 'cell', False): (nonlinear_jacobian_apply_volume_templates, nonlinear_jacobian_apply_volume_args), - ('jacobian_apply', 'exterior_facet', False): (nonlinear_jacobian_apply_boundary_templates, nonlinear_jacobian_apply_boundary_args), - ('jacobian_apply', 'interior_facet', False): (nonlinear_jacobian_apply_skeleton_templates, nonlinear_jacobian_apply_skeleton_args), - }.get((form_type, integral_type, linear), None) + ('jacobian_apply', 'exterior_facet', True): (jacobian_apply_boundary_templates, jacobian_apply_boundary_args), + ('jacobian_apply', 'interior_facet', True): (jacobian_apply_skeleton_templates, jacobian_apply_skeleton_args), + ('jacobian_apply', 'cell', False): (nonlinear_jacobian_apply_volume_templates, nonlinear_jacobian_apply_volume_args), + ('jacobian_apply', 'exterior_facet', False): (nonlinear_jacobian_apply_boundary_templates, nonlinear_jacobian_apply_boundary_args), + ('jacobian_apply', 'interior_facet', False): (nonlinear_jacobian_apply_skeleton_templates, nonlinear_jacobian_apply_skeleton_args), + }.get((form_type, integral_type, linear), None) return construct_signature(templates(), args(), kernel_name()) +def assembly_routine_args(): + integral_type = get_global_context_value("integral_type") + form_type = get_global_context_value("form_type") + formdata = get_global_context_value("formdata") + + args = {('residual', 'cell'): alpha_volume_args, + ('residual', 'exterior_facet'): alpha_boundary_args, + ('residual', 'interior_facet'): alpha_skeleton_args, + ('jacobian', 'cell'): jacobian_volume_args, + ('jacobian', 'exterior_facet'): jacobian_boundary_args, + ('jacobian', 'interior_facet'): jacobian_skeleton_args, + }.get((form_type, integral_type), None) + + if args is None: + # Check if form is linear + from dune.perftool.pdelab.driver import is_linear + linear = is_linear(formdata.original_form) + + args = {('jacobian_apply', 'cell', True): jacobian_apply_volume_args, + ('jacobian_apply', 'exterior_facet', True): jacobian_apply_boundary_args, + ('jacobian_apply', 'interior_facet', True): jacobian_apply_skeleton_args, + ('jacobian_apply', 'cell', False): nonlinear_jacobian_apply_volume_args, + ('jacobian_apply', 'exterior_facet', False): nonlinear_jacobian_apply_boundary_args, + ('jacobian_apply', 'interior_facet', False): nonlinear_jacobian_apply_skeleton_args, + }.get((form_type, integral_type, linear), None) + + return args() + + def construct_signature(types, args, name): templates = "template<{}>".format(", ".join("typename {}".format(t) for t in set(types))) func = "void {}({}) const".format(name, ", ".join("{}{}& {}".format("const " if c else "", t, a) for t, (c, a) in zip(types, args))) diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py index 1bbcf245888c29526bb46a7de6c848110e40e1bb..d82ea08af1a7ad439562be9732114b098848d6c1 100644 --- a/python/dune/perftool/sumfact/__init__.py +++ b/python/dune/perftool/sumfact/__init__.py @@ -8,6 +8,7 @@ from dune.perftool.sumfact.basis import (lfs_inames, pymbolic_trialfunction, pymbolic_trialfunction_gradient, ) +import dune.perftool.sumfact.switch from dune.perftool.pdelab import PDELabInterface diff --git a/python/dune/perftool/sumfact/switch.py b/python/dune/perftool/sumfact/switch.py new file mode 100644 index 0000000000000000000000000000000000000000..78b31590f8ea3656b19254ec77dde955c4607982 --- /dev/null +++ b/python/dune/perftool/sumfact/switch.py @@ -0,0 +1,107 @@ +""" boundary and skeleton integrals come in variants in sum factorization - implement the switch! """ + +from dune.perftool.generation import (backend, + get_global_context_value, + global_context, + ) +from dune.perftool.pdelab.geometry import world_dimension +from dune.perftool.pdelab.localoperator import generate_kernel +from dune.perftool.pdelab.signatures import (assembly_routine_args, + assembly_routine_signature, + kernel_name, + ) +from dune.perftool.cgen.clazz import ClassMember + + +@backend(interface="generate_kernels_per_integral", name="sumfact") +def generate_kernels_per_integral(integrals): + dim = get_global_context_value("formdata").geometric_dimension + measure = get_global_context_value("integral_type") + + if measure == "cell": + yield generate_kernel(integrals) + + if measure == "exterior_facet": + # Generate all necessary kernels + for facedir in range(dim): + for facemod in range(2): + with global_context(facedir_s=facedir, facemod_s=facemod): + yield generate_kernel(integrals) + + # Generate switch statement + yield generate_exterior_facet_switch() + + if measure == "interior_facet": + # Generate all necessary kernels + for facedir_s in range(dim): + for facemod_s in range(2): + for facedir_n in range(dim): + for facemod_n in range(2): + with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n): + yield generate_kernel(integrals) + + # Generate switch statement + yield generate_interior_facet_switch() + + +def get_kernel_name(facedir_s=None, facemod_s=None, facedir_n=None, facemod_n=None): + with global_context(facedir_s=facedir_s, facemod_s=facemod_s, facedir_n=facedir_n, facemod_n=facemod_n): + return kernel_name() + + +def generate_exterior_facet_switch(): + # Extract the signature + signature = assembly_routine_signature() + args = ", ".join(tuple(a for c, a in assembly_routine_args())) + dim = world_dimension() + + # Construct the switch statement + block = [] + block.append("{") + block.append(" size_t variant = ig.indexInInside();") + block.append(" switch(variant)") + block.append(" {") + + for facedir_s in range(dim): + for facemod_s in range(2): + block.append(" case {}: {}({}); break;".format(dim * facedir_s + facemod_s, + get_kernel_name(facedir_s=facedir_s, + facemod_s=facemod_s, + ), + args)) + + block.append(" }") + block.append("}") + + return ClassMember(signature + block) + + +def generate_interior_facet_switch(): + # Extract the signature + signature = assembly_routine_signature() + args = ", ".join(tuple(a for c, a in assembly_routine_args())) + dim = world_dimension() + + # Construct the switch statement + block = [] + block.append("{") + block.append(" size_t variant = ig.indexInOutside() + 6 * ig.indexInInside();") + block.append(" switch(variant)") + block.append(" {") + + for facedir_s in range(dim): + for facemod_s in range(2): + for facedir_n in range(dim): + for facemod_n in range(2): + block.append(" case {}: {}({}); break;".format((dim * facedir_s + facemod_s) * (2 * dim) + dim * facedir_n + facemod_n, + get_kernel_name(facedir_s=facedir_s, + facemod_s=facemod_s, + facedir_n=facedir_n, + facemod_n=facemod_n, + ), + args)) + + block.append(" }") + block.append("}") + + return ClassMember(signature + block)