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.

# 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
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!

(Please note I do not further address the correctness of your variational form)
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  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »