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   
My approach reads:

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
***
***     fenics-support@googlegroups.com
***
*** 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  

1 Answer


0
4 months ago by
It seems that periodic BCs and DG are still incompatible in the current version of FEniCS

https://bitbucket.org/fenics-project/dolfin/issues/558/periodic-boundary-conditions-for-dg
Please login to add an answer/comment or follow this question.

Similar posts:
Search »