### How to correctly implement integral over facets in DG?

393

views

1

I am trying to implement the follwing linear variation problem:

$\int_{\Omega_e}\phi^{n+1}vdx=\int_{\Omega_e}\phi^nvdx-\Delta t\int_{\Omega_e}\nabla\cdot\left(u\phi^n\right)vdx+\Delta t\int_{\partial\Omega_e}\frac{1}{2}\left(n\cdot u-\left|n\cdot u\right|\right)\left(\phi^n-\phi^{n,+}\right)vdS$∫_{Ωe}`ϕ`^{n+1}`v``d``x`=∫_{Ωe}`ϕ`^{n}`v``d``x`−Δ`t`∫_{Ωe}∇·(`u``ϕ`^{n})`v``d``x`+Δ`t`∫_{∂Ωe}12 (`n`·`u`−|`n`·`u`|)(`ϕ`^{n}−`ϕ`^{n,+})`v``d``S` My implementation, inspired by https://github.com/FEniCS/dolfin/blob/master/demo/undocumented/dg-advection-diffusion/python/demo_dg-advection-diffusion.py:

```
from dolfin import *
# FIXME: Make mesh ghosted
parameters["ghost_mode"] = "shared_facet"
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
# return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)
return bool((near(x[0], 0) or near(x[1], 0)) and
(not ((near(x[0], 0) and near(x[1], 1)) or
(near(x[0], 1) and near(x[1], 0)))) and on_boundary)
def map(self, x, y):
if near(x[0], 1) and near(x[1], 1):
y[0] = x[0] - 1.
y[1] = x[1] - 1.
elif near(x[0], 1):
y[0] = x[0] - 1.
y[1] = x[1]
else: # near(x[1], 1)
y[0] = x[0]
y[1] = x[1] - 1.
# Load mesh
mesh = UnitSquareMesh( 6, 6 )
# Defining the function spaces
V_dg = FunctionSpace(mesh, "DG", 4)
V_cg = FunctionSpace(mesh, "CG", 4)
V_u = VectorFunctionSpace(mesh, "CG", 2)
# Create velocity Function
uForm = Constant( ( "1.0", "1.0" ) )
u = project( uForm, V_u )
# Test and trial functions
v = TestFunction(V_dg)
phi = TrialFunction(V_dg)
# Mesh-related functions
n = FacetNormal(mesh)
un = (dot(u, n) - abs(dot(u, n)))/2.0
initial = Expression( "sin( 2 * pi * x[ 0 ] ) * sin( 2 * pi * x[ 1 ] )", element = V_dg.ufl_element() )
phi0 = project( initial, V_dg )
dt = 0.0001
t = 0
while t < 0.1:
t = t + dt
# Bilinear form
a = phi * v * dx \
- dt * dot( u * phi, grad( v ) ) * dx \
+ dt * un('-') * ( phi0('-') - phi0('+') ) * v('-') * dS
# Linear form
L = phi0 * v * dx
# Solution function
phi_h = Function(V_dg)
# solve the system
solve( a == L, phi_h )
# swap generations
phi0 = phi_h
# Project solution to a continuous function space
up = project(phi_h, V=V_cg)
file = File("temperature.pvd")
file << up
```

This results in the error:

```
raise ArityMismatch("Integrand arguments {0} differ from form arguments {1}.".format(args, arguments))
ArityMismatch: Integrand arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Discontinuous Lagrange', triangle, 4)), 0, None),) differ from form arguments (Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Discontinuous Lagrange', triangle, 4)), 0, None), Argument(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', triangle, 1), dim=2), 0), FiniteElement('Discontinuous Lagrange', triangle, 4)), 1, None)).
```

It only occurs if the integral over the facets are included. How can I solve this problem?

Community: FEniCS Project

### 1 Answer

4

Hi,

you get this error for the simple reason that the facet term

So to avoid running into this error, you can modify your code as

Hope this helps!

(Please note I do not further address the correctness of your variational form)

you get this error for the simple reason that the facet term

`+ dt * un('-') * ( phi0('-') - phi0('+') ) * v('-') * dS`

factually is a linear form. So to avoid running into this error, you can modify your code as

```
# Bilinear form
a = phi * v * dx \
- dt * dot( u * phi, grad( v ) ) * dx # \
# Linear form
L = phi0 * v * dx \
- dt * un('-') * ( phi0('-') - phi0('+') ) * v('-') * dS
```

Hope this helps!

(Please note I do not further address the correctness of your variational form)

2

I am indeed not quite confident about your facet integral term. A thesis which I always find useful in developing FEniCS DG formulations is the work by Oelgaard (in particular Chapter 4).

written
8 months ago by
Jakob Maljaars

Thank you for the help and reference!

written
8 months ago by
Philipp O

Please login to add an answer/comment or follow this question.

yes it does. Thank you!

Did you find a problem with the variational form? I am happy for input. Especially as my implementation appears to be instable.