diff --git a/python/dune/perftool/blockstructured/geometry.py b/python/dune/perftool/blockstructured/geometry.py index 14d0b6f6785051b5767898b1ed05e521638eecf7..0fe5e00c37ae75e978e0a17f5baed437c6560ba9 100644 --- a/python/dune/perftool/blockstructured/geometry.py +++ b/python/dune/perftool/blockstructured/geometry.py @@ -1,16 +1,15 @@ from dune.perftool.generation import (get_backend, temporary_variable, - instruction, - preamble) + instruction) from dune.perftool.tools import get_pymbolic_basename from dune.perftool.options import (get_form_option, option_switch) -from dune.perftool.pdelab.geometry import (name_jacobian_determinant, - world_dimension, +from dune.perftool.pdelab.geometry import (world_dimension, local_dimension, apply_to_global_transformation, name_facet_jacobian_determinant, - name_element_corners) + name_element_corners, + component_iname) from dune.perftool.blockstructured.tools import sub_element_inames import pymbolic.primitives as prim from loopy.match import Writes @@ -18,14 +17,49 @@ from loopy.match import Writes # TODO warum muss within_inames_is_final=True gesetzt werden? def compute_jacobian_2d(name): - pos = get_backend("quad_pos", selector=option_switch("blockstructured"))() + pymbolic_pos = get_backend("quad_pos", selector=option_switch("blockstructured"))() corners = name_element_corners() - instruction(expression=prim.Subscript(prim.Variable(corners), (0,0)), - assignee=prim.Subscript(prim.Variable(name), (0,0)), - within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()), - within_inames_is_final=True, - depends_on=frozenset({Writes(corners), Writes(get_pymbolic_basename(pos))}) - ) + + jac_inames = tuple(component_iname(context="jac", count=i) for i in range(world_dimension())) + + def _corner_diff(first, second, additional_terms=tuple()): + simplified_names = tuple("cd" + n.split('_')[2] + n.split('_')[3] for n in additional_terms) + name = "corner_diff_"+"_".join((str(first), str(second))+simplified_names) + temporary_variable(name, shape_impl=('fv',), shape=(world_dimension(),)) + diminame = component_iname(context='corner') + + if additional_terms: + xs_sum = prim.Sum(tuple(prim.Subscript(prim.Variable(x), (prim.Variable(diminame),)) for x in additional_terms)) + else: + xs_sum = 0 + + instruction(expression=prim.Sum((prim.Subscript(prim.Variable(corners), (first, prim.Variable(diminame))), + -1 * prim.Subscript(prim.Variable(corners), (second, prim.Variable(diminame))), + -1 * xs_sum)), + assignee=prim.Subscript(prim.Variable(name), prim.Variable(diminame)), + within_inames=frozenset(), within_inames_is_final=True, + depends_on=frozenset({Writes(corners)})) + return name + + c = _corner_diff(2, 0) + b = _corner_diff(1, 0) + a = _corner_diff(3, 0, (b, c)) + + expr_jac = [None, None] + expr_jac[0] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (1,)), + prim.Subscript(prim.Variable(a), (prim.Variable(jac_inames[0]),)))), + prim.Subscript(prim.Variable(b), (prim.Variable(jac_inames[0]),)))) + expr_jac[1] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (0,)), + prim.Subscript(prim.Variable(a), (prim.Variable(jac_inames[1]),)))), + prim.Subscript(prim.Variable(c), (prim.Variable(jac_inames[1]),)))) + + for i, expr in enumerate(expr_jac): + instruction(expression=expr, + assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_inames[i]),i)), + within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()), + within_inames_is_final=True, + depends_on=frozenset({Writes(b), Writes(get_pymbolic_basename(pymbolic_pos))}) + ) @@ -44,14 +78,18 @@ def name_jacobian_matrix(): def pymbolic_compute_determinant(matrix): - pass + dim = world_dimension() + matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)] + expr_determinant = prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])), + -1*prim.Product((matrix_entry[1][0], matrix_entry[0][1])))) + return expr_determinant def define_jacobian_determinant(name): temporary_variable(name, shape=()) jacobian = name_jacobian_matrix() - determinant = pymbolic_compute_determinant(jacobian) - return instruction(expression=prim.Subscript(prim.Variable(jacobian), (0,0)), + expr_determinant = pymbolic_compute_determinant(jacobian) + return instruction(expression=prim.Call(prim.Variable("abs"), (expr_determinant,)), assignee=prim.Variable(name), within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()), within_inames_is_final=True, @@ -67,7 +105,13 @@ def name_jacobian_determinant(): # scale determinant according to the order of the blockstructure def pymbolic_jacobian_determinant(): - return prim.Quotient(prim.Variable(name_jacobian_determinant()), + if get_form_option("constant_transformation_matrix"): + from dune.perftool.pdelab.geometry import name_jacobian_determinant + n_det = name_jacobian_determinant() + else: + from dune.perftool.blockstructured.geometry import name_jacobian_determinant + n_det = name_jacobian_determinant() + return prim.Quotient(prim.Variable(n_det), prim.Power(get_form_option("number_of_blocks"), local_dimension()))