From e7470548f16153e4df148e9cba3388dfa0926df5 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Tue, 31 Jan 2017 16:54:24 +0100
Subject: [PATCH] Implement handlers to allow generation of correct penalty
term
---
.../convection_diffusion/conv_diff_dg.ufl | 14 ++++++---
applications/poisson_dg/poisson_dg.ufl | 14 ++++++---
python/dune/perftool/loopy/target.py | 2 +-
python/dune/perftool/pdelab/__init__.py | 8 +++++
python/dune/perftool/pdelab/geometry.py | 30 +++++++++++++++++++
python/dune/perftool/ufl/preprocess.py | 18 +++++------
python/dune/perftool/ufl/visitor.py | 16 ++++++++++
test/poisson/opcount_poisson_dg.ufl | 12 +++++---
test/poisson/poisson_dg.ufl | 16 ++++++----
test/poisson/poisson_dg_neumann.ufl | 16 ++++++----
test/poisson/poisson_dg_tensor.ufl | 16 ++++++----
11 files changed, 124 insertions(+), 38 deletions(-)
diff --git a/applications/convection_diffusion/conv_diff_dg.ufl b/applications/convection_diffusion/conv_diff_dg.ufl
index 38e8dc73..80e92055 100644
--- a/applications/convection_diffusion/conv_diff_dg.ufl
+++ b/applications/convection_diffusion/conv_diff_dg.ufl
@@ -1,3 +1,4 @@
+dim = 3
x = SpatialCoordinate(cell)
g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
@@ -13,17 +14,22 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
theta = -1.0
r = (inner(A*grad(u), grad(v)) + (c*u-f)*v)*dx \
+ inner(n, A*avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
+ theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
- inner(n, A*grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
- theta*u*inner(A*grad(v), n)*ds \
+ theta*g*inner(A*grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/applications/poisson_dg/poisson_dg.ufl b/applications/poisson_dg/poisson_dg.ufl
index 53636053..7d532811 100644
--- a/applications/poisson_dg/poisson_dg.ufl
+++ b/applications/poisson_dg/poisson_dg.ufl
@@ -1,3 +1,4 @@
+dim = 3
x = SpatialCoordinate(cell)
f = -6.
g = x[0]*x[0] + x[1]*x[1] + x[2]*x[2]
@@ -9,18 +10,23 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
+
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
+ theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
- theta*u*inner(grad(v), n)*ds \
- f*v*dx \
+ theta*g*inner(grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/python/dune/perftool/loopy/target.py b/python/dune/perftool/loopy/target.py
index 9427a401..6ef27648 100644
--- a/python/dune/perftool/loopy/target.py
+++ b/python/dune/perftool/loopy/target.py
@@ -118,7 +118,7 @@ class DuneCExpressionToCodeMapper(CExpressionToCodeMapper):
result = self.rec(children.pop(), PREC_NONE)
while children:
- result = "std::%s(%s, %s)" % (what,
+ result = "%s(%s, %s)" % (what,
self.rec(children.pop(), PREC_NONE),
result)
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index 8529eb1e..dcc17068 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -11,6 +11,8 @@ from dune.perftool.pdelab.basis import (pymbolic_basis,
pymbolic_reference_gradient,
)
from dune.perftool.pdelab.geometry import (dimension_iname,
+ name_cell_volume,
+ name_facet_area,
name_facet_jacobian_determinant,
name_jacobian_determinant,
name_jacobian_inverse_transposed,
@@ -112,6 +114,12 @@ class PDELabInterface(object):
def name_unit_outer_normal(self):
return name_unit_outer_normal()
+ def name_cell_volume(self, restriction):
+ return name_cell_volume(restriction)
+
+ def name_facet_area(self):
+ return name_facet_area()
+
#
# Quadrature related generator functions
#
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index b66f58fc..ce080eba 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -405,3 +405,33 @@ def to_global(local):
name = get_pymbolic_basename(local) + "_global"
apply_to_global_transformation(name, local)
return prim.Variable(name)
+
+
+def define_cell_volume(name, restriction):
+ geo = name_cell_geometry(restriction)
+ temporary_variable(name, shape=())
+ code = "{} = {}.volume();".format(name, geo)
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ )
+
+
+def name_cell_volume(restriction):
+ name = restricted_name("volume", restriction)
+ define_cell_volume(name, restriction)
+ return name
+
+
+def define_facet_area(name):
+ geo = name_intersection_geometry()
+ temporary_variable(name, shape=())
+ code = "{} = {}.volume();".format(name, geo)
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ )
+
+
+def name_facet_area():
+ name = "area"
+ define_facet_area(name)
+ return name
diff --git a/python/dune/perftool/ufl/preprocess.py b/python/dune/perftool/ufl/preprocess.py
index 540300bb..913bf5da 100644
--- a/python/dune/perftool/ufl/preprocess.py
+++ b/python/dune/perftool/ufl/preprocess.py
@@ -1,22 +1,20 @@
""" Preprocessing algorithms for UFL forms """
+import ufl.classes as uc
-def preprocess_form(form):
- from ufl.classes import (FacetJacobianDeterminant,
- FacetNormal,
- JacobianDeterminant,
- JacobianInverse,
- )
+def preprocess_form(form):
from ufl.algorithms import compute_form_data
formdata = compute_form_data(form,
do_apply_function_pullbacks=True,
do_apply_geometry_lowering=True,
do_apply_integral_scaling=True,
- preserve_geometry_types=(FacetJacobianDeterminant,
- FacetNormal,
- JacobianDeterminant,
- JacobianInverse,
+ preserve_geometry_types=(uc.CellVolume,
+ uc.FacetArea,
+ uc.FacetJacobianDeterminant,
+ uc.FacetNormal,
+ uc.JacobianDeterminant,
+ uc.JacobianInverse,
)
)
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index 1ba817c8..de92f3a4 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -265,6 +265,12 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
raise NotImplementedError("Power function not really implemented")
+ def max_value(self, o):
+ return prim.Max(tuple(self.call(op) for op in o.ufl_operands))
+
+ def min_value(self, o):
+ return prim.Min(tuple(self.call(op) for op in o.ufl_operands))
+
#
# Handler for conditionals, use pymbolic base implementation
#
@@ -360,6 +366,16 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
def facet_jacobian_determinant(self, o):
return Variable(self.interface.name_facet_jacobian_determinant())
+ def cell_volume(self, o):
+ restriction = self.restriction
+ if self.measure == 'exterior_facet':
+ restriction = Restriction.NEGATIVE
+
+ return prim.Variable(self.interface.name_cell_volume(restriction))
+
+ def facet_area(self, o):
+ return prim.Variable(self.interface.name_facet_area())
+
#
# Equality/Hashability of the visitor
# In order to pass this visitor into (cached) generator functions, it is beneficial
diff --git a/test/poisson/opcount_poisson_dg.ufl b/test/poisson/opcount_poisson_dg.ufl
index 3ced1b82..ed537bc9 100644
--- a/test/poisson/opcount_poisson_dg.ufl
+++ b/test/poisson/opcount_poisson_dg.ufl
@@ -11,20 +11,24 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
+ theta*u*inner(grad(v), n)*ds \
- f*v*dx \
- theta*g*inner(grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/test/poisson/poisson_dg.ufl b/test/poisson/poisson_dg.ufl
index daf2793e..b8ad0425 100644
--- a/test/poisson/poisson_dg.ufl
+++ b/test/poisson/poisson_dg.ufl
@@ -1,11 +1,13 @@
cell = triangle
+degree = 1
+dim = 2
x = SpatialCoordinate(cell)
c = (0.5-x[0])**2 + (0.5-x[1])**2
g = exp(-1.*c)
f = 4*(1.-c)*g
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
@@ -13,20 +15,24 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
+ theta*u*inner(grad(v), n)*ds \
- f*v*dx \
- theta*g*inner(grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
diff --git a/test/poisson/poisson_dg_neumann.ufl b/test/poisson/poisson_dg_neumann.ufl
index 1dd15d4f..884787e6 100644
--- a/test/poisson/poisson_dg_neumann.ufl
+++ b/test/poisson/poisson_dg_neumann.ufl
@@ -1,3 +1,5 @@
+dim = 3
+degree = 1
cell = triangle
x = SpatialCoordinate(cell)
@@ -8,7 +10,7 @@ sgn = conditional(x[1] > 0.5, 1., -1.)
j = -2.*sgn*(x[1]-0.5)*g
bctype = conditional(Or(x[1]<1e-8, x[1]>1.-1e-8), 0, 1)
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
@@ -18,21 +20,25 @@ n = FacetNormal(cell)('+')
ds = ds(subdomain_data=bctype)
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds(1) \
- + gamma*u*v*ds(1) \
+ + gamma_ext*u*v*ds(1) \
+ theta*u*inner(grad(v), n)*ds(1) \
- f*v*dx \
- theta*g*inner(grad(v), n)*ds(1) \
- - gamma*g*v*ds(1) \
+ - gamma_ext*g*v*ds(1) \
- j*v*ds(0)
forms = [r]
diff --git a/test/poisson/poisson_dg_tensor.ufl b/test/poisson/poisson_dg_tensor.ufl
index d84cebf6..1ae0139b 100644
--- a/test/poisson/poisson_dg_tensor.ufl
+++ b/test/poisson/poisson_dg_tensor.ufl
@@ -1,3 +1,5 @@
+dim = 2
+degree = 1
cell = quadrilateral
x = SpatialCoordinate(cell)
@@ -8,7 +10,7 @@ g = x[0]**2 + x[1]**2
c = 8.
f = -4.
-V = FiniteElement("DG", cell, 1)
+V = FiniteElement("DG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
@@ -16,20 +18,24 @@ v = TestFunction(V)
n = FacetNormal(cell)('+')
# penalty factor
-gamma = 1.0
+alpha = 1.0
+h_ext = CellVolume(cell) / FacetArea(cell)
+gamma_ext = (alpha * degree * (degree + dim - 1)) / h_ext
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = (inner(A*grad(u), grad(v)) + c*u*v)*dx \
+ inner(n, A*avg(grad(u)))*jump(v)*dS \
- + gamma*jump(u)*jump(v)*dS \
+ + gamma_int*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(A*avg(grad(v)), n)*dS \
- inner(n, A*grad(u))*v*ds \
- + gamma*u*v*ds \
+ + gamma_ext*u*v*ds \
+ theta*u*inner(A*grad(v), n)*ds \
- f*v*dx \
- theta*g*inner(A*grad(v), n)*ds \
- - gamma*g*v*ds
+ - gamma_ext*g*v*ds
forms = [r]
--
GitLab