diff --git a/python/dune/perftool/blockstructured/geometry.py b/python/dune/perftool/blockstructured/geometry.py
index c43ac2c18b0ef425d00ba28f1c3b6ab2f092b426..627f09daae7a96922a265731a3e9e66b3a6b4d75 100644
--- a/python/dune/perftool/blockstructured/geometry.py
+++ b/python/dune/perftool/blockstructured/geometry.py
@@ -17,50 +17,90 @@ import pymbolic.primitives as prim
from loopy.match import Writes
-# TODO warum muss within_inames_is_final=True gesetzt werden?
-def compute_jacobian_2d(name):
- pymbolic_pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+def compute_corner_diff(first, second, additional_terms=tuple()):
corners = name_element_corners()
+ 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')
- jac_inames = tuple(component_iname(context="jac", count=i) for i in range(world_dimension()))
+ 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)} | {Writes(term) for term in additional_terms}))
+ return name
- 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
+# TODO warum muss within_inames_is_final=True gesetzt werden?
+# Transformation T(x,y) = a * xy + b * x + c * y + d
+def compute_jacobian_2d(name):
+ pymbolic_pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+ jac_iname = tuple(component_iname(context="jac", count=i) for i in range(world_dimension()))
- c = _corner_diff(2, 0)
- b = _corner_diff(1, 0)
- a = _corner_diff(3, 0, (b, c))
+ c = compute_corner_diff(2, 0)
+ b = compute_corner_diff(1, 0)
+ a = compute_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]),))))
+ prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)))),
+ prim.Subscript(prim.Variable(b), (prim.Variable(jac_iname),))))
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]),))))
+ prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)))),
+ prim.Subscript(prim.Variable(c), (prim.Variable(jac_iname),))))
+
+ for i, expr in enumerate(expr_jac):
+ instruction(expression=expr,
+ assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_iname),i)),
+ within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()),
+ within_inames_is_final=True,
+ depends_on=frozenset({Writes(a), Writes(b), Writes(c), Writes(get_pymbolic_basename(pymbolic_pos))})
+ )
+
+
+# Transformation T(x,y,z) = a * xyz + b * xy + c * xz + d * yz + e * x + f * y + g * z + h
+def compute_jacobian_3d(name):
+ pymbolic_pos = get_backend("quad_pos", selector=option_switch("blockstructured"))()
+ jac_iname = component_iname(context="jac")
+
+ g = compute_corner_diff(4, 0)
+ f = compute_corner_diff(2, 0)
+ e = compute_corner_diff(1, 0)
+ d = compute_corner_diff(6, 0, (f, g))
+ c = compute_corner_diff(5, 0, (e, g))
+ b = compute_corner_diff(3, 0, (e, f))
+ a = compute_corner_diff(7, 0, (b, c, d, e, f, g))
+ all_cds = [a, b, c, d, e, f, g]
+
+ expr_jac = [None, None, None]
+ terms = [[b, c, e], [b, d, f], [c, d, g]]
+
+ # j[:][i] = a * qp[k]*qp[l] + v1 * qp[k] + v2 * qp[l] + v3
+ # with k, l in {0,1,2} != i and k<l and vj = terms[i][j]
+ for i in range(3):
+ k, l = sorted(set(range(3))-{i})
+ expr_jac[i] = prim.Sum((prim.Product((prim.Subscript(pymbolic_pos, (k,)), prim.Subscript(pymbolic_pos, (l,)),
+ prim.Subscript(prim.Variable(a), (prim.Variable(jac_iname),)))),
+ prim.Product((prim.Subscript(pymbolic_pos, (k,)),
+ prim.Subscript(prim.Variable(terms[i][0]), (prim.Variable(jac_iname),)))),
+ prim.Product((prim.Subscript(pymbolic_pos, (l,)),
+ prim.Subscript(prim.Variable(terms[i][1]), (prim.Variable(jac_iname),)))),
+ prim.Subscript(prim.Variable(terms[i][2]), (prim.Variable(jac_iname),))
+ ))
for i, expr in enumerate(expr_jac):
instruction(expression=expr,
- assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_inames[i]),i)),
+ assignee=prim.Subscript(prim.Variable(name), (prim.Variable(jac_iname),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))})
+ depends_on=frozenset({Writes(get_pymbolic_basename(pymbolic_pos))}).union(frozenset([Writes(cd) for cd in all_cds]))
)
@@ -68,6 +108,8 @@ def define_jacobian_matrix(name):
temporary_variable(name, shape_impl=('fm',), shape=(world_dimension(), world_dimension()))
if world_dimension() == 2:
compute_jacobian_2d(name)
+ elif world_dimension() == 3:
+ compute_jacobian_3d(name)
else:
raise NotImplementedError()
@@ -80,18 +122,26 @@ def name_jacobian_matrix():
def compute_determinant(name, matrix):
dim = world_dimension()
+ matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
if dim == 2:
- 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]))))
- instruction(expression=expr_determinant,
- assignee=prim.Variable(name),
- within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()),
- within_inames_is_final=True,
- depends_on=frozenset({Writes(matrix)})
- )
+ elif dim == 3:
+ expr_determinant = prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1], matrix_entry[2][2])),
+ prim.Product((matrix_entry[0][1], matrix_entry[1][2], matrix_entry[2][0])),
+ prim.Product((matrix_entry[0][2], matrix_entry[1][0], matrix_entry[2][1])),
+ -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][1], matrix_entry[2][0])),
+ -1 * prim.Product((matrix_entry[0][0], matrix_entry[1][2], matrix_entry[2][1])),
+ -1 * prim.Product((matrix_entry[0][1], matrix_entry[1][0], matrix_entry[2][2]))
+ ))
else:
raise NotImplementedError()
+ instruction(expression=expr_determinant,
+ assignee=prim.Variable(name),
+ within_inames=frozenset(sub_element_inames()+get_backend(interface="quad_inames")()),
+ within_inames_is_final=True,
+ depends_on=frozenset({Writes(matrix)})
+ )
def define_jacobian_determinant(name):