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

393
views
1
8 months ago by

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+1vdx=ΩeϕnvdxΔtΩe·(uϕn)vdx+Δt∂Ωe12 (n·u|n·u|)(ϕnϕn,+)vdS 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.

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

4
8 months ago by
Hi,

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!

Hi,

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.
written 8 months ago by Philipp O
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