Error or bug in implementation of DG for 2D advection equation

278
views
0
8 months ago by

I am trying to implement the following DG-formulation of a 2D advection equation with periodic boundary conditions:
$\int_{\Omega_e}\left(\frac{\partial}{\partial t}f+\nabla\cdot\left(cf\right)\right)\phi dx=\int_{\partial\Omega_e}\frac{1}{2}\left(n\cdot c-\left|n\cdot c\right|\right)\left(f-f^+\right)\phi ds$Ωe(t ƒ +·(cƒ ))ϕdx=∂Ωe12 (n·c|n·c|)(ƒ ƒ +)ϕds
Here,  $n$n denotes the outer normal and  $c$c  denotes a constant vector.
For testing, I want to use explicit Euler for time discretization:
$\int_{\Omega_e}f^{n+1}\phi dx=\int_{\Omega_e}f^n\phi dx-\Delta t\int_{\Omega_e}\nabla\cdot\left(cf^n\right)\phi dx+\Delta t\int_{\partial\Omega_e}\frac{1}{2}\left(n\cdot c-\left|n\cdot c\right|\right)\left(f-f^+\right)\phi ds$Ωeƒ n+1ϕdx=Ωeƒ nϕdxΔtΩe·(cƒ n)ϕdx+Δt∂Ωe12 (n·c|n·c|)(ƒ ƒ +)ϕds

from dolfin import *

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.

dt = 0.5

mesh = UnitSquareMesh( 20, 20 )

pbc = PeriodicBoundary()
Scalar = FunctionSpace( mesh, "DG", 3, constrained_domain = pbc )

v = TestFunction( Scalar )
u = TrialFunction( Scalar )

e_alpha = Constant( ( 1, 1 ) )
n = FacetNormal( mesh )

initial = Expression( "sin( 2 * x[ 0 ] * pi ) * sin( 2 * x[ 1 ] * pi )", element = Scalar.ufl_element() )
u0 = project( initial, Scalar )

m = u * v * dx # mass matrix for time t_{n+1}
rhs = u0 * v * dx \ # mass matrix for time t_n
- dt * dot( e_alpha, grad( u0 ) ) * v * dx \ # advection term
+ dt * 0.5 * ( dot( n, e_alpha ) - abs( dot( n, e_alpha ) ) ) * jump( u0 ) * v * ds # integral over edges

ures = Function( Scalar )
solve( m == rhs, ures )

This crashes with error message:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to perform just-in-time compilation of form.
*** Reason:  ffc.jit failed with message:
Traceback (most recent call last):
File "/Users/philipp/anaconda/envs/fenicsproject/lib/python2.7/site-packages/dolfin/compilemodules/jit.py", line 142, in jit
result = ffc.jit(ufl_object, parameters=p)
File "/Users/philipp/anaconda/envs/fenicsproject/lib/python2.7/site-packages/ffc/jitcompiler.py", line 210, in jit
raise FFCJitError("A directory with files to reproduce the jit build failure has been created.")

The compilation using recompile.sh results in error message:

ffc_form_f01bab02420d23eaca8e9473ec308dd691c1e135.cpp:282:51: error: use of undeclared identifier 'facet_1'; did you mean 'facet'?
w1_r1 += w[1][ic + 10] * FE10_C0_F_Q4[facet_1][iq][ic];
^~~~~~~
facet
ffc_form_f01bab02420d23eaca8e9473ec308dd691c1e135.cpp:193:49: note: 'facet' declared here
std::size_t facet,
^
ffc_form_f01bab02420d23eaca8e9473ec308dd691c1e135.cpp:285:46: error: use of undeclared identifier 'facet_0'; did you mean 'facet'?
w1_r0 += w[1][ic] * FE10_C0_F_Q4[facet_0][iq][ic];
^~~~~~~
facet
ffc_form_f01bab02420d23eaca8e9473ec308dd691c1e135.cpp:193:49: note: 'facet' declared here
std::size_t facet,
^
2 errors generated.


Did I implement the variational form incorrectly or is this a bug?

Community: FEniCS Project
Aren't you missing terms on the interior facets dS?

Edit: Also UFL/FEniCS sadly does not recognise that jumps and averages take different values depending on whether you're on an interior or exterior facet. You'll have to write these each explicitly.

There's a DG advection-diffusion demo. I recommend you examine that for tips.
written 8 months ago by Nate
Thank you for pointing to dS. Using dS instead of dS already helped a lot.
written 8 months ago by Philipp O