From 66ef042c323225e1032a7f1f64b0c9780dfd7aaa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Mon, 6 Aug 2018 15:56:32 +0200
Subject: [PATCH] Sum factorized evaluation of geometry mappings
Calculate the evaluation and jacobians of geometry mappings of cells
using sum factorization. By avoiding all geometry mapping of intersections
we can avoid the face-twist problem and just make the grid consistens (not
yet done).
TODO:
- Make grid consistent
- Vectorization
- 3D
---
dune/perftool/sumfact/invertgeometry.hh | 19 +
.../dune/perftool/blockstructured/__init__.py | 6 +-
.../dune/perftool/blockstructured/geometry.py | 2 +-
python/dune/perftool/pdelab/__init__.py | 6 +-
python/dune/perftool/pdelab/geometry.py | 61 +--
python/dune/perftool/pdelab/quadrature.py | 1 -
python/dune/perftool/sumfact/__init__.py | 34 +-
python/dune/perftool/sumfact/accumulation.py | 1 -
python/dune/perftool/sumfact/geometry.py | 472 +++++++++++++++---
python/dune/perftool/sumfact/quadrature.py | 86 +---
python/dune/perftool/sumfact/symbolic.py | 17 +
python/dune/perftool/ufl/visitor.py | 2 +-
12 files changed, 505 insertions(+), 202 deletions(-)
create mode 100644 dune/perftool/sumfact/invertgeometry.hh
diff --git a/dune/perftool/sumfact/invertgeometry.hh b/dune/perftool/sumfact/invertgeometry.hh
new file mode 100644
index 00000000..a09f5e52
--- /dev/null
+++ b/dune/perftool/sumfact/invertgeometry.hh
@@ -0,0 +1,19 @@
+#ifndef DUNE_PERFTOOL_SUMFACT_INVERTGEOMETRY_HH
+#define DUNE_PERFTOOL_SUMFACT_INVERTGEOMETRY_HH
+
+
+template<typename T>
+inline T invert_and_return_determinant(const T a00, const T a10, const T a01, const T a11, T jacobian_inverse[4]){
+ T det = a00 * a11 - a10 * a01;
+ assert (det != 0.0);
+
+ jacobian_inverse[0] = a11 / det;
+ jacobian_inverse[1] = -a10 / det;
+ jacobian_inverse[2] = -a01 / det;
+ jacobian_inverse[3] = a00 / det;
+
+ return std::abs(det);
+}
+
+
+#endif
diff --git a/python/dune/perftool/blockstructured/__init__.py b/python/dune/perftool/blockstructured/__init__.py
index 8849491e..5c1283b2 100644
--- a/python/dune/perftool/blockstructured/__init__.py
+++ b/python/dune/perftool/blockstructured/__init__.py
@@ -9,7 +9,7 @@ from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position
from dune.perftool.blockstructured.spaces import lfs_inames
from dune.perftool.blockstructured.basis import (pymbolic_reference_gradient,
pymbolic_basis)
-from dune.perftool.blockstructured.geometry import (pymbolic_jacobian_inverse_transposed,
+from dune.perftool.blockstructured.geometry import (pymbolic_jacobian_inverse,
pymbolic_jacobian_determinant,
pymbolic_facet_jacobian_determinant,
to_global)
@@ -55,8 +55,8 @@ class BlockStructuredInterface(PDELabInterface):
def pymbolic_jacobian_determinant(self):
return pymbolic_jacobian_determinant()
- def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
- return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ def pymbolic_jacobian_inverse(self, i, j, restriction):
+ return pymbolic_jacobian_inverse(i, j, restriction)
def pymbolic_facet_jacobian_determinant(self):
return pymbolic_facet_jacobian_determinant()
diff --git a/python/dune/perftool/blockstructured/geometry.py b/python/dune/perftool/blockstructured/geometry.py
index 1ffbdf04..113f5940 100644
--- a/python/dune/perftool/blockstructured/geometry.py
+++ b/python/dune/perftool/blockstructured/geometry.py
@@ -19,7 +19,7 @@ def pymbolic_jacobian_determinant():
# scale Jacobian according to the order of the blockstructure
-def pymbolic_jacobian_inverse_transposed(i, j, restriction):
+def pymbolic_jacobian_inverse(i, j, restriction):
name_jit = get_backend(interface="name_jit", selector=option_switch("constant_transformation_matrix"))(restriction)
return prim.Product((get_form_option("number_of_blocks"),
prim.Subscript(prim.Variable(name_jit), (j, i))))
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index d3153ca0..f53b6ff6 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -19,7 +19,7 @@ from dune.perftool.pdelab.geometry import (component_iname,
pymbolic_facet_area,
pymbolic_facet_jacobian_determinant,
pymbolic_jacobian_determinant,
- pymbolic_jacobian_inverse_transposed,
+ pymbolic_jacobian_inverse,
pymbolic_unit_inner_normal,
pymbolic_unit_outer_normal,
to_global,
@@ -125,8 +125,8 @@ class PDELabInterface(object):
def pymbolic_jacobian_determinant(self):
return pymbolic_jacobian_determinant()
- def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
- return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ def pymbolic_jacobian_inverse(self, i, j, restriction):
+ return pymbolic_jacobian_inverse(i, j, restriction)
def pymbolic_unit_inner_normal(self):
return pymbolic_unit_inner_normal()
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 39a46f21..16f75495 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -204,14 +204,14 @@ def name_in_cell_geometry(restriction):
# TODO is it always necessary to add quadrature inames?
def apply_in_cell_transformation(name, local, restriction):
geo = name_in_cell_geometry(restriction)
- return get_backend("quadrature_preamble")("{} = {}.global({});".format(name,
- geo,
- str(local),
- ),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(local)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(local))}),
- )
+ return quadrature_preamble("{} = {}.global({});".format(name,
+ geo,
+ str(local),
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(local)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(local))}),
+ )
def pymbolic_in_cell_coordinates(local, restriction):
@@ -253,11 +253,11 @@ def local_dimension():
def evaluate_unit_outer_normal(name):
ig = name_intersection_geometry_wrapper()
qp = get_backend("quad_pos")()
- return get_backend("quadrature_preamble")("{} = {}.unitOuterNormal({});".format(name, ig, qp),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(qp)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
- )
+ return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
+ )
@preamble
@@ -267,7 +267,7 @@ def declare_normal(name, kernel, decl_info):
def pymbolic_unit_outer_normal():
- name = "outer_normal"
+ name = "unit_outer_normal"
if not get_form_option("diagonal_transformation_matrix"):
temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
evaluate_unit_outer_normal(name)
@@ -340,14 +340,14 @@ def define_jacobian_inverse_transposed(name, restriction):
geo = name_cell_geometry(restriction)
pos = get_backend("qp_in_cell", selector=option_switch(("blockstructured", "sumfact")))(restriction)
- return get_backend("quadrature_preamble")("{} = {}.jacobianInverseTransposed({});".format(name,
- geo,
- str(pos),
- ),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
+ return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
+ geo,
+ str(pos),
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
@backend(interface="name_jit", name="default")
@@ -357,10 +357,11 @@ def name_jacobian_inverse_transposed(restriction):
return name
-def pymbolic_jacobian_inverse_transposed(i, j, restriction):
+def pymbolic_jacobian_inverse(i, j, restriction):
+ # Calculate jacobian_inverse_transposed
name = get_backend(interface="name_jit", selector=option_switch("constant_transformation_matrix"))(restriction)
- # Dune only has JacobianInverseTransposed as a first class citizen,
- # so we need to switch the indices around!
+
+ # Return jacobian inverse -> (j, i) instead of (i, j)
return prim.Subscript(prim.Variable(name), (j, i))
@@ -397,11 +398,11 @@ def define_jacobian_determinant(name):
str(pos),
)
- return get_backend("quadrature_preamble")(code,
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
@backend(interface="fdetjac", name="constant_transformation_matrix")
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 61008799..031d97b0 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -61,7 +61,6 @@ def quadrature_inames():
return (quadrature_iname(),)
-@backend(interface="quadrature_preamble")
def quadrature_preamble(code, **kw):
kw['tags'] = kw.get('tags', frozenset({})).union(frozenset({"quad"}))
return instruction(inames=get_backend(interface="quad_inames")(), code=code, **kw)
diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py
index d5eed554..d3192b6f 100644
--- a/python/dune/perftool/sumfact/__init__.py
+++ b/python/dune/perftool/sumfact/__init__.py
@@ -78,37 +78,25 @@ class SumFactInterface(PDELabInterface):
return ret
def pymbolic_unit_outer_normal(self):
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
- ret, indices = pymbolic_unit_outer_normal(self.visitor.indices)
- self.visitor.indices = indices
- return ret
+ from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
+ return pymbolic_unit_outer_normal(self.visitor)
def pymbolic_unit_inner_normal(self):
from dune.perftool.sumfact.geometry import pymbolic_unit_inner_normal
- ret, indices = pymbolic_unit_inner_normal(self.visitor.indices)
- self.visitor.indices = indices
- return ret
+ return pymbolic_unit_inner_normal(self.visitor)
def pymbolic_facet_jacobian_determinant(self):
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
- return pymbolic_facet_jacobian_determinant()
+ from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
+ return pymbolic_facet_jacobian_determinant(self.visitor)
def pymbolic_facet_area(self):
from dune.perftool.sumfact.geometry import pymbolic_facet_area
- return pymbolic_facet_area()
+ return pymbolic_facet_area(self.visitor)
- def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
- from dune.perftool.pdelab.geometry import pymbolic_jacobian_inverse_transposed
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ def pymbolic_jacobian_inverse(self, i, j, restriction):
+ from dune.perftool.sumfact.geometry import pymbolic_jacobian_inverse
+ return pymbolic_jacobian_inverse(i, j, restriction, self.visitor)
def pymbolic_jacobian_determinant(self):
- from dune.perftool.pdelab.geometry import pymbolic_jacobian_determinant
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- return pymbolic_jacobian_determinant()
+ from dune.perftool.sumfact.geometry import pymbolic_jacobian_determinant
+ return pymbolic_jacobian_determinant(self.visitor)
diff --git a/python/dune/perftool/sumfact/accumulation.py b/python/dune/perftool/sumfact/accumulation.py
index 4e7c2a69..abad6788 100644
--- a/python/dune/perftool/sumfact/accumulation.py
+++ b/python/dune/perftool/sumfact/accumulation.py
@@ -24,7 +24,6 @@ from dune.perftool.options import (get_form_option,
)
from dune.perftool.loopy.flatten import flatten_index
from dune.perftool.loopy.target import type_floatingpoint
-from dune.perftool.sumfact.quadrature import nest_quadrature_loops
from dune.perftool.pdelab.driver import FEM_name_mangling
from dune.perftool.pdelab.localoperator import determine_accumulation_space
from dune.perftool.pdelab.restriction import restricted_name
diff --git a/python/dune/perftool/sumfact/geometry.py b/python/dune/perftool/sumfact/geometry.py
index f3e5b9fd..ecec31af 100644
--- a/python/dune/perftool/sumfact/geometry.py
+++ b/python/dune/perftool/sumfact/geometry.py
@@ -6,38 +6,71 @@ from dune.perftool.generation import (backend,
get_backend,
get_counted_variable,
iname,
+ include_file,
instruction,
kernel_cached,
preamble,
+ silenced_warning,
temporary_variable,
+ get_global_context_value,
globalarg,
+ valuearg,
)
+from dune.perftool.options import get_option
from dune.perftool.pdelab.geometry import (local_dimension,
world_dimension,
+ name_cell_geometry,
name_geometry,
)
-from dune.perftool.sumfact.switch import get_facedir
-from dune.perftool.sumfact.symbolic import SumfactKernelInterfaceBase
+from dune.perftool.pdelab.localoperator import (name_ansatz_gfs_constructor_param,
+ lop_template_ansatz_gfs,
+ lop_template_range_field,
+ )
+from dune.perftool.pdelab.restriction import restricted_name
+from dune.perftool.sumfact.basis import construct_basis_matrix_sequence
+from dune.perftool.sumfact.quadrature import (additional_inames,
+ default_quadrature_inames)
+from dune.perftool.sumfact.realization import (name_buffer_storage,
+ realize_sum_factorization_kernel,
+ )
+from dune.perftool.sumfact.switch import get_facedir, get_facemod
+from dune.perftool.sumfact.symbolic import SumfactKernelInterfaceBase, SumfactKernel
+from dune.perftool.sumfact.tabulation import (quadrature_points_per_direction,
+ BasisTabulationMatrix,
+ )
from dune.perftool.sumfact.vectorization import attach_vectorization_info
+from dune.perftool.tools import get_pymbolic_basename
from dune.perftool.options import get_form_option, option_switch
from dune.perftool.ufl.modified_terminals import Restriction
from pytools import ImmutableRecord
+from loopy.match import Writes
+
import pymbolic.primitives as prim
import numpy as np
@iname
-def corner_iname():
- name = get_counted_variable("corneriname")
- domain(name, 2 ** local_dimension())
+def global_corner_iname(restriction):
+ name = get_counted_variable(restricted_name("global_corneriname", restriction))
+ domain(name, 2 ** world_dimension())
return name
class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
- def __init__(self, dir):
- ImmutableRecord.__init__(self, dir=dir)
+ def __init__(self, direction, restriction):
+ """Base class for sum-factorized evaluation of geometry mappings
+
+ At the moment we only do this for cells and not faces. For
+ intersections we do this corresponding reference elements of the
+ neigboring cells.
+
+ Each spatial component needs a seperate sum factorization kernel. The
+ argument 'direction' specifies the component (x-component: 0,
+ y-component: 1, z-component: 2).
+ """
+ ImmutableRecord.__init__(self, direction=direction, restriction=restriction)
@property
def stage(self):
@@ -48,18 +81,22 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
def realize(self, sf, insn_dep):
# TODO: When we use vectorization this should be the offset
+ assert(get_form_option('vectorization_strategy') == 'none')
index = 0
- from dune.perftool.sumfact.realization import name_buffer_storage
+ # Note: world_dimension, since we only do evaluation of cell geometry mappings
name = "input_{}".format(sf.buffer)
temporary_variable(name,
- shape=(2 ** local_dimension(), sf.vector_width),
+ shape=(2 ** world_dimension(), sf.vector_width),
custom_base_storage=name_buffer_storage(sf.buffer, 0),
managed=True,
)
+ ciname = global_corner_iname(self.restriction)
- ciname = corner_iname()
- geo = name_geometry()
+ if self.restriction == 0:
+ geo = name_geometry()
+ else:
+ geo = name_cell_geometry(self.restriction)
# NB: We need to realize this as a C instruction, because the corner
# method does return a non-scalar, which does not fit into the current
@@ -71,7 +108,7 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
index,
geo,
ciname,
- self.dir,
+ self.direction,
)
insn = instruction(code=code,
@@ -85,23 +122,33 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
@backend(interface="spatial_coordinate", name="default")
def pymbolic_spatial_coordinate_multilinear(do_predicates, visitor):
+ """Calculate global coordinates of quadrature points for multilinear geometry mappings
+
+ The evalualation is done using sum factorization techniques.
+
+ On facets we use the geometry mapping of the self cell to get a global
+ evaluation of the quadrature points. Avoiding the geometry mapping of the
+ face itself greatly simplifies the generation of code for unstructured
+ grids. Instead of looking at all the possible embeddings of the reference
+ element of the face into the reference element of the cell we can make the
+ grid edge consistent and choose a suitable convention for the order of
+ directions in our sum factorization kernels.
+ """
assert len(visitor.indices) == 1
- # Construct the matrix sequence for the evaluation of the global coordinate.
- # We need to manually construct this one, because on facets, we want to use the
- # geometry embedding of the facet into the global space directly without going
- # through the neighboring cell geometries. That matrix sequence will only have
- # dim-1 matrices!
- from dune.perftool.sumfact.tabulation import quadrature_points_per_direction, BasisTabulationMatrix
- quadrature_size = quadrature_points_per_direction()
- matrix_sequence = tuple(BasisTabulationMatrix(direction=i, basis_size=2) for i in range(world_dimension()) if i != get_facedir(visitor.restriction))
- inp = GeoCornersInput(visitor.indices[0])
-
- from dune.perftool.sumfact.symbolic import SumfactKernel
+ # Adjust restriction for boundary facets
+ restriction = visitor.restriction
+ if get_global_context_value("integral_type") == "exterior_facet":
+ restriction = 1
+
+ # Generate sum factorization kernel and add vectorization info
+ matrix_sequence = construct_basis_matrix_sequence(facedir=get_facedir(restriction),
+ facemod=get_facemod(restriction),
+ basis_size=(2,) * world_dimension())
+ inp = GeoCornersInput(visitor.indices[0], restriction)
sf = SumfactKernel(matrix_sequence=matrix_sequence,
interface=inp,
)
-
vsf = attach_vectorization_info(sf)
# Add a sum factorization kernel that implements the evaluation of
@@ -109,7 +156,10 @@ def pymbolic_spatial_coordinate_multilinear(do_predicates, visitor):
from dune.perftool.sumfact.realization import realize_sum_factorization_kernel
var, _ = realize_sum_factorization_kernel(vsf)
+ # The results of this function is already the right component of the
+ # spatial coordinate and we don't need to index this in the visitor
visitor.indices = None
+
return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
@@ -164,7 +214,6 @@ def pymbolic_spatial_coordinate_axiparallel(do_predicates, visitor):
# Urgh: *SOMEHOW* construct a face direction
from dune.perftool.pdelab.restriction import Restriction
restriction = Restriction.NONE
- from dune.perftool.generation import get_global_context_value
if get_global_context_value("integral_type") == "interior_facet":
restriction = Restriction.POSITIVE
from dune.perftool.sumfact.switch import get_facedir
@@ -192,69 +241,171 @@ def pymbolic_spatial_coordinate_axiparallel(do_predicates, visitor):
return prim.Subscript(prim.Variable(lowcorner), (index,)) + x * prim.Subscript(prim.Variable(meshwidth), (index,))
-def pymbolic_unit_outer_normal(visitor_indices):
- index, = visitor_indices
+def define_outer_normal(name, visitor):
+ """Calculate outer normal (length not necessarily 1)
+
+ This will be used to calculate the unit outer normal and the determinant of
+ the Jacobian of the geometry mapping of the face.
+
+ The outer normal is calculated by multiplying the Jacobian inverse
+ transposed of the geometry mapping of the 'self'-cell with the unit outer
+ normal \hat{n_s} of the face corresponding to the intersetcion in the
+ reference element of the 'self'-cell:
+
+ n = (\nabla \mu_{T_s})^{-T} \hat{n_s}
+
+ Instead of setting up \hat{n_s} we just look at facedir_s and facemod_s.
+ """
+ facedir_s = get_facedir(Restriction.POSITIVE)
+ facemod_s = get_facemod(Restriction.POSITIVE)
+ temporary_variable(name, shape=(world_dimension(),))
+ for i in range(world_dimension()):
+ assignee = prim.Subscript(prim.Variable(name), i)
+ jit = pymbolic_jacobian_inverse(i, facedir_s, Restriction.POSITIVE, visitor)
+
+ # Note: 2*facemod_s-1 because of
+ # -1 if facemod_s = 0
+ # +1 if facemod_s = 1
+ expression = jit * (2 * facemod_s - 1)
+
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ )
+
+
+def name_outer_normal(visitor):
+ name = "outer_normal"
+ define_outer_normal(name, visitor)
+ return name
+
+
+def pymbolic_outer_normal(visitor):
+ name = name_outer_normal(visitor)
+ return prim.Variable(name)
+
+
+def define_norm_of_outer_normal(name, visitor):
+ outer_normal = pymbolic_outer_normal(visitor)
+
+ # Calculate the norm of outer_normal
+ temporary_variable(name, shape=())
+ expression = 0
+ for i in range(world_dimension()):
+ expression = prim.Sum((expression, prim.Subscript(outer_normal, i) * prim.Subscript(outer_normal, i)))
+ expression = prim.Call(prim.Variable("sqrt"), (expression,))
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+
+ # Apparently the call of sqrt shadows the dependency on outer normal. Add manually.
+ instruction(assignee=prim.Variable(name),
+ expression=expression,
+ depends_on=frozenset({Writes(get_pymbolic_basename(outer_normal))}),
+ within_inames=frozenset(inames),
+ )
+
+
+def name_norm_of_outer_normal(visitor):
+ name = "norm_of_outer_normal"
+ define_norm_of_outer_normal(name, visitor)
+ return name
+
+
+def pymbolic_norm_of_outer_normal(visitor):
+ name = name_norm_of_outer_normal(visitor)
+ return prim.Variable(name)
+
+
+def pymbolic_unit_outer_normal(visitor):
+ assert (len(visitor.indices) is 1)
+ index = visitor.indices[0]
assert isinstance(index, int)
if get_form_option("diagonal_transformation_matrix"):
+ # In this case we just return -1, 0 or 1 depending on the current
+ # facemod and facedir. We don't need to (and cannot) index these ints
+ # so we set the indices of the visitor to None.
+ visitor.indices = None
+
+ # Use facemod_s and facedir_s
from dune.perftool.sumfact.switch import get_facedir, get_facemod
if index == get_facedir(Restriction.POSITIVE):
if get_facemod(Restriction.POSITIVE):
- return 1, None
+ return 1
else:
- return -1, None
+ return -1
else:
- return 0, None
+ return 0
else:
- from dune.perftool.pdelab.geometry import pymbolic_unit_outer_normal as _norm
- return _norm(), visitor_indices
+ # In pymbolic_outer_normal we need to create the
+ # jacobian_inverse_transposed which doesn't play nicely with the
+ # visitor having indices. Maybe this should be improved. For now let's
+ # just set them to None and restore them afterwards.
+ visitor.indices = None
+ outer_normal = pymbolic_outer_normal(visitor)
+ norm = pymbolic_norm_of_outer_normal(visitor)
+ visitor.indices = (index,)
+
+ # Create unit_outer_normal by scaling outer_normal to length 1
+ name = "unit_outer_normal"
+ temporary_variable(name, shape=(world_dimension(),))
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ for i in range(world_dimension()):
+ instruction(assignee=prim.Subscript(prim.Variable(name), i),
+ expression=prim.Subscript(outer_normal, i) / norm,
+ within_inames=frozenset(inames),
+ )
+
+ return prim.Variable(name)
def pymbolic_unit_inner_normal(visitor_indices):
- index, = visitor_indices
+ # Note: The implementation of this was adjusted for unstructured grids but
+ # it was not properly tested. If you need the unit inner normal just remove
+ # the assert(False) statement but make sure to verify that this does the
+ # right thing.
+ assert(False)
+
+ assert (len(visitor.indices) is 1)
+ index = visitor.indices[0]
assert isinstance(index, int)
if get_form_option("diagonal_transformation_matrix"):
- from dune.perftool.sumfact.switch import get_facedir, get_facemod
+ # In this case we just return -1, 0 or 1 depending on the current
+ # facemod and facedir. We don't need to (and cannot) index these ints
+ # so we set the indices of the visitor to None.
+ visitor.indices = None
if index == get_facedir(Restriction.POSITIVE):
if get_facemod(Restriction.POSITIVE):
- return -1, None
+ return -1
else:
- return 1, None
+ return 1
else:
- return 0, None
- else:
- from dune.perftool.pdelab.geometry import pymbolic_unit_inner_normal as _norm
- return _norm(), visitor_indices
-
-
-def pymbolic_facet_jacobian_determinant():
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_constant_facet_jacobian_determinant()
+ return 0
else:
- from dune.perftool.pdelab.geometry import pymbolic_facet_jacobian_determinant as _norm
- return _norm()
+ # In pymbolic_outer_normal we need to create the
+ # jacobian_inverse_transposed which doesn't play nicely with the
+ # visitor having indices. Maybe this should be improved. For now let's
+ # just set them to None and restore them afterwards.
+ visitor.indices = None
+ outer_normal = pymbolic_outer_normal(visitor)
+ norm = pymbolic_norm_of_outer_normal(visitor)
+ visitor.indices = (index,)
+ # Create unit_outer_normal by scaling outer_normal to length 1
+ name = "unit_outer_normal"
+ temporary_variable(name, shape=(world_dimension(),))
-def pymbolic_constant_facet_jacobian_determinant():
- facedir = get_facedir(Restriction.POSITIVE)
- assert isinstance(facedir, int)
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ for i in range(world_dimension()):
+ instruction(assignee=prim.Subscript(prim.Variable(name), i),
+ expression=-1 * prim.Subscript(outer_normal, i) / norm,
+ within_inames=frozenset(inames),
+ )
- name = "fdetjac"
- define_constant_facet_jacobian_determinant(name)
- globalarg(name, shape=(world_dimension(),))
- return prim.Subscript(prim.Variable(name), (facedir,))
-
-
-@class_member(classtag="operator")
-def define_constant_facet_jacobian_determinant(name):
- define_constant_facet_jacobian_determinant_eval(name)
- from dune.perftool.pdelab.localoperator import lop_template_ansatz_gfs
- gfst = lop_template_ansatz_gfs()
- return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
+ return prim.Variable(name)
@kernel_cached(kernel="operator")
def define_constant_facet_jacobian_determinant_eval(name):
- from dune.perftool.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
gfs = name_ansatz_gfs_constructor_param()
rft = lop_template_range_field()
code = ["{",
@@ -269,9 +420,196 @@ def define_constant_facet_jacobian_determinant_eval(name):
)
-def pymbolic_facet_area():
+@class_member(classtag="operator")
+def define_constant_facet_jacobian_determinant(name):
+ define_constant_facet_jacobian_determinant_eval(name)
+ gfst = lop_template_ansatz_gfs()
+ return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
+
+
+def pymbolic_constant_facet_jacobian_determinant():
+ facedir = get_facedir(Restriction.POSITIVE)
+ assert isinstance(facedir, int)
+
+ name = "fdetjac"
+ define_constant_facet_jacobian_determinant(name)
+ globalarg(name, shape=(world_dimension(),))
+ return prim.Subscript(prim.Variable(name), (facedir,))
+
+
+def define_facet_jacobian_determinant(name, visitor):
+ """Absolute value of determinant of jacobian of facet geometry mapping
+
+ Calculate the absolute value of the determinant of the jacobian of the
+ geometry mapping from the reference cell of the intersection to the
+ intersection:
+
+ |\det(\nabla \mu_F)|
+
+ This is done using the surface metric tensor lemma from the lecture notes
+ of Jean-Luc Guermond:
+
+ |\det(\nabla \mu_F)| = \| \tilde{n} \| |\det(\nabla \mu_{T_s})|
+
+ Here \tilde{n} is the outer normal defined by:
+
+ \tilde{n} = (\nabla \mu_{T_s})^{-T} \hat{n}_s
+
+ where \hat{n}_s is the unit outer normal of the corresponding face of the
+ reference cell.
+ """
+ # Calculate norm of outer normal and determinant of jacobian of geometry
+ # mapping of the inside cell -> set restriction to 1 for boundary facets
+ if get_global_context_value("integral_type") == "exterior_facet":
+ assert visitor.restriction is 0
+ visitor.restriction = 1
+ norm_of_outer_normal = pymbolic_norm_of_outer_normal(visitor)
+ detjac = pymbolic_jacobian_determinant(visitor)
+ if get_global_context_value("integral_type") == "exterior_facet":
+ visitor.restriction = 0
+
+ # Multiply detjac and norm_of_outer_normal
+ temporary_variable(name, shape=())
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ instruction(assignee=prim.Variable(name),
+ expression=detjac * norm_of_outer_normal,
+ within_inames=frozenset(inames),
+ )
+
+
+def name_facet_jacobian_determinant(visitor):
+ name = "fdetjac"
+ define_facet_jacobian_determinant(name, visitor)
+ return name
+
+
+def pymbolic_facet_jacobian_determinant(visitor):
+ if get_form_option("constant_transformation_matrix"):
+ return pymbolic_constant_facet_jacobian_determinant()
+ else:
+ name = name_facet_jacobian_determinant(visitor)
+ return prim.Variable(name)
+
+
+def pymbolic_facet_area(visitor):
+ if get_form_option("constant_transformation_matrix"):
+ return pymbolic_facet_jacobian_determinant(visitor)
+ else:
+ from dune.perftool.pdelab.geometry import pymbolic_facet_area as nonconstant_facet_area
+ return nonconstant_facet_area()
+
+
+def _name_jacobian(i, j, restriction, visitor):
+ """Return the (i, j) component of the jacobian of the geometry mapping
+
+ Evaluation of the derivative of the geometry mapping is done using sum
+ factorization.
+
+ Note: At the moment this only works for the mappings from reference cells
+ to the cell and not for the geometry mappings of intersections.
+ """
+ do_predicates = visitor.do_predicates
+
+ # Create matrix sequence with derivative in j direction
+ matrix_sequence = construct_basis_matrix_sequence(derivative=j,
+ facedir=get_facedir(restriction),
+ facemod=get_facemod(restriction),
+ basis_size=(2,) * world_dimension())
+
+ # Sum factorization input for the i'th component of the geometry mapping
+ inp = GeoCornersInput(i, restriction)
+
+ sf = SumfactKernel(matrix_sequence=matrix_sequence,
+ interface=inp,
+ )
+ vsf = attach_vectorization_info(sf)
+
+ # Add a sum factorization kernel that implements the evaluation of
+ # the basis functions at quadrature points (stage 1)
+ var, _ = realize_sum_factorization_kernel(vsf)
+
+ assert(visitor.indices is None)
+ return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
+
+
+def define_jacobian_inverse(name, restriction, visitor):
+ """Return jacobian inverse transposed of the geometry mapping (of a cell)
+
+ At the moment this only works for geometry mappings of cells and not for
+ intersection. We only consider this case as it greatly simplifies code
+ generation for unstructured grids by making the grid edge consistent and
+ avoiding the face-twist problem.
+ """
+ # Calculate the jacobian of the geometry mapping using sum factorization
+ dim = world_dimension()
+ names_jacobian = []
+ for j in range(dim):
+ for i in range(dim):
+ name_jacobian = restricted_name("jacobian_geometry_{}{}".format(i, j), restriction)
+ temporary_variable(name_jacobian)
+ assignee = prim.Variable(name_jacobian)
+ expression = _name_jacobian(i, j, restriction, visitor)
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ )
+ names_jacobian.append(name_jacobian)
+
+ # Calculate the inverse of the jacobian of the geometry mapping and the
+ # determinant by calling a c++ function
+ name_detjac = name_jacobian_determinant(visitor)
+ temporary_variable(name_detjac, shape=())
+ temporary_variable(name, shape=(dim, dim), managed=True)
+ include_file('dune/perftool/sumfact/invertgeometry.hh', filetag='operatorfile')
+ code = "{} = invert_and_return_determinant({}, {});".format(name_detjac,
+ ", ".join(names_jacobian),
+ name)
+ silenced_warning("read_no_write({})".format(name_detjac))
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ return instruction(code=code,
+ assignees=frozenset([name, name_detjac]),
+ read_variables=frozenset(names_jacobian),
+ inames=frozenset(inames),
+ )
+
+
+def name_jacobian_inverse(restriction, visitor):
+ name = restricted_name("jacobian_inverse", restriction)
+ define_jacobian_inverse(name, restriction, visitor)
+ return name
+
+
+def pymbolic_jacobian_inverse(i, j, restriction, visitor):
+ # Switch between implementations without using the backend mechanism. In
+ # the case of constant_transformation_matrix we want to use the pdelab
+ # branch, otherwise we want to go into the sumfact implementation instead
+ # of pdelab.
if get_form_option("constant_transformation_matrix"):
- return pymbolic_facet_jacobian_determinant()
+ from dune.perftool.pdelab.geometry import name_constant_jacobian_inverse_transposed
+ name = name_constant_jacobian_inverse_transposed(restriction)
+ # Return jacobian inverse -> (j, i) instead of (i, j)
+ return prim.Subscript(prim.Variable(name), (j, i))
else:
- from dune.perftool.pdelab.geometry import pymbolic_facet_area as _norm
- return _norm()
+ name = name_jacobian_inverse(restriction, visitor)
+ return prim.Subscript(prim.Variable(name), (i, j))
+
+
+def name_jacobian_determinant(visitor):
+ name = "detjac"
+ return name
+
+
+def pymbolic_jacobian_determinant(visitor):
+ if get_form_option("constant_transformation_matrix"):
+ from dune.perftool.pdelab.geometry import pymbolic_jacobian_determinant
+ return pymbolic_jacobian_determinant()
+ else:
+ # The calculation of the jacobian determinant happens in this function
+ if visitor.measure == 'cell':
+ restriction = 0
+ else:
+ restriction = 1
+ name_jacobian_inverse(restriction, visitor)
+
+ return prim.Variable(name_jacobian_determinant(visitor))
diff --git a/python/dune/perftool/sumfact/quadrature.py b/python/dune/perftool/sumfact/quadrature.py
index 586d9fd2..6fe70fa4 100644
--- a/python/dune/perftool/sumfact/quadrature.py
+++ b/python/dune/perftool/sumfact/quadrature.py
@@ -39,20 +39,6 @@ import loopy as lp
import numpy as np
-def nest_quadrature_loops(kernel, inames):
- insns = []
- for insn in kernel.instructions:
- if "quad" in insn.tags:
- insns.append(insn.copy(within_inames=insn.within_inames.union(frozenset(inames)),
- tags=insn.tags - frozenset({"quad"})
- )
- )
- else:
- insns.append(insn)
-
- return kernel.copy(instructions=insns)
-
-
class BaseWeight(FunctionIdentifier):
def __init__(self, accumobj):
self.accumobj = accumobj
@@ -80,6 +66,15 @@ def pymbolic_base_weight():
return 1.0
+def default_quadrature_inames(visitor):
+ info = visitor.current_info[1]
+ if info is None:
+ element = None
+ else:
+ element = info.element
+ return quadrature_inames(element)
+
+
@kernel_cached
def quadrature_inames(element):
if element is None:
@@ -177,7 +172,8 @@ def quadrature_weight(visitor):
element = None
else:
element = info.element
- return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
+ return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"),
+ tuple(prim.Variable(i) for i in quadrature_inames(element)))
def define_quadrature_position(name, local_index):
@@ -245,7 +241,9 @@ def pymbolic_indexed_quadrature_position(local_index, visitor):
return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
-def _additional_inames(visitor):
+def additional_inames(visitor):
+ """Return inames for iterating over ansatz space in the case of jacobians
+ """
info = visitor.current_info[1]
if info is None:
element = None
@@ -262,59 +260,3 @@ def _additional_inames(visitor):
return lfs_inames
else:
return ()
-
-
-@backend(interface="quadrature_preamble", name="sumfact")
-def quadrature_preamble(code, **kw):
- visitor = get_global_context_value('visitor')
-
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
-
- inames = quadrature_inames(element)
- inames = inames + _additional_inames(visitor)
-
- return instruction(inames=inames, code=code, **kw)
-
-
-@preamble
-def define_quadrature_point(name):
- local_dim = local_dimension()
- return "Dune::FieldVector<double, {}> {};".format(local_dim, name)
-
-
-def name_quadrature_point():
- name = "quadrature_point"
- define_quadrature_point(name)
- return name
-
-
-@backend(interface="quad_pos", name="sumfact")
-def pymbolic_quadrature_position():
- visitor = get_global_context_value('visitor')
-
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
-
- quad_point = name_quadrature_point()
- inames = quadrature_inames(element) + _additional_inames(visitor)
- globalarg(quad_point, shape=(local_dimension(),))
- facedir = get_facedir(visitor.restriction)
- for i in range(local_dimension()):
- expression = pymbolic_indexed_quadrature_position(i, visitor)
- instruction(assignee=prim.Subscript(prim.Variable(quad_point), i),
- expression=expression,
- within_inames=frozenset(inames))
- return Variable(quad_point)
-
-
-@backend(interface="qp_in_cell", name="sumfact")
-def pymbolic_quadrature_position_in_cell(restriction):
- from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
- return pymbolic_quadrature_position_in_cell(restriction)
diff --git a/python/dune/perftool/sumfact/symbolic.py b/python/dune/perftool/sumfact/symbolic.py
index fb283a05..f2fbcb76 100644
--- a/python/dune/perftool/sumfact/symbolic.py
+++ b/python/dune/perftool/sumfact/symbolic.py
@@ -417,6 +417,23 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
if len(self.matrix_sequence) == local_dimension():
return tuple(prim.Variable(i) for i in quad_inames)
+ # # Extremly ugly. This is the special case of being on an intersection
+ # # and generating a SF kernel for the geometry evaluation of one of the
+ # # neighboring cells. This screws the usual setup since the SF sequence
+ # # is not equal to the local dimension.
+ # #
+ # # TODO: Vectorization is probably broken for this case!
+ # count_sf_sequence_length = 0
+ # for d in range(world_dimension()):
+ # if self.matrix_sequence[d].face is not None:
+ # count_sf_sequence_length += 1
+ # if count_sf_sequence_length != local_dimension():
+ # from dune.perftool.sumfact.quadrature import global_quadrature_inames
+ # quad_inames = global_quadrature_inames(element)
+
+ # for d in range(world_dimension()):
+ # print(self.matrix_sequence[d].face)
+
# Traverse all the quadrature inames and map them to their correct direction
index = []
i = 0
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index cbede6c3..9e14391c 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -464,7 +464,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
if isinstance(i, int) and isinstance(j, int) and i != j:
return 0
- return self.interface.pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ return self.interface.pymbolic_jacobian_inverse(i, j, restriction)
def jacobian(self, o):
raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")
--
GitLab