# DG variational form

**URL:** <https://fenicsproject.discourse.group/t/dg-variational-form/742>\
**Category:** Uncategorized\
**Created:** [May 6, 2019, 6:43am UTC](https://fenicsproject.discourse.group/t/dg-variational-form/742 "2019-05-06T06:43:03Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![reddeviln](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/reddeviln/32/164_2.png) [@reddeviln](https://fenicsproject.discourse.group/u/reddeviln)\
**Post date:** [May 6, 2019, 6:43am UTC](https://fenicsproject.discourse.group/t/dg-variational-form/742/1 "2019-05-06T06:43:03Z")

</div>

Hi everyone,

the problem I am trying to solve works fine when using FEM but I am having trouble rewriting it to use it with DG elements. Maybe I have something wrong about the fluxes. The following are the equations I am trying to solve:  
U\_t=\begin{pmatrix}u\\C\_1\\C\_2\end{pmatrix}\_t = F\_x(U) +G\_y(U)+q(U)  
where F=\begin{pmatrix}C\_1 -1\\\frac{1}{\lambda}u\\0\end{pmatrix}, G=\begin{pmatrix}C\_2 -1 \\0\\\frac{1}{\lambda}u\end{pmatrix} are the fluxes in the 2 dimensional domain and q=\begin{pmatrix}-\nabla p\\\frac{-1}{\lambda}(C\_1-1)\\\frac{-1}{\lambda}(C\_1-1)\end{pmatrix} the source term. Where \lambda is a coefficient and \nabla p is a given pressure gradient. Now when I derive the variational form with a testfunction \phi i get  
\int U\_t\phi d\Omega -\int F\phi n\_xd\partial\Omega-\int G\phi n\_yd\partial\Omega-\int F\phi\_xd\Omega-\int G\phi\_yd\Omega=\int q\phi d\Omega where n is the cell normal. I use the local lax friedrichs flux for the numerical flux terms: f(u\_+,u\_-)=\frac{1}{2}(f(u\_+)+f(u\_-)) - \frac{1}{2\sqrt{\lambda}}(u\_+-u\_-) where I have put the maximum Eigenvalue in already. The solver for my problem runs just fine. But it is unstable however it should be stable. My code:

```python
from FEniCSSimulation import *
import mshr
import matplotlib.pyplot as plt

#init and mesh
nCells = 40
domain=mshr.Circle(Point(0,0),sqrt(9))
mesh = mshr.generate_mesh(domain, nCells)
WithStressTensorSim.mesh = mesh
# dofs
StressTensorElement1 = FiniteElement('DG', triangle, 0)
StressTensorElement2 = FiniteElement('DG', triangle, 0)
VelocityElements = FiniteElement('DG', triangle, 0)
element = MixedElement(StressTensorElement1, StressTensorElement2, VelocityElements)
V = FunctionSpace(mesh, element)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V.sub(2), Constant(0), boundary)
u_n = interpolate(Constant((0,0,0)), V)
Lambda = 1
dt= Constant(0.1)
n =FacetNormal(self.mesh)
U =TrialFunction(self.V[0])
C1, C2, u = split(self.U)
testf = TestFunction(self.V[0])
CV1, CV2, v = split(self.testf)
un1, un2 , un3 = split(u_n)
#numerical flux terms beyond here
Fnum1int = avg(C1-1)
Fnum1ext = (C1 - 1) - 1 / (2 * sqrt(Lambda)) * u
Fnum2int = 1 / Lambda * avg(u)
Fnum2ext = 1 / Lambda * u - 1 / (2 * sqrt(Lambda)) * C1
Fnum3int = 0
Fnum3ext = - 1 / (2 * sqrt(Lambda)) * C2
Gnum1int = avg(C2 - 1)
Gnum1ext = (C2 - 1) - 1 / (2 * sqrt(Lambda)) * u
Gnum2int = 0
Gnum2ext = - 1 / (2 * sqrt(Lambda)) * C1
Gnum3int = 1 / Lambda * avg(u)
Gnum3ext = 1 / Lambda * u - 1 / (2 * sqrt(Lambda)) * C2

result = dot(as_vector([Fnum1int, Gnum1int]) - 1 / (2 * sqrt(Lambda)) * jump(u, n), jump(v, n)) * dS \
            + dot(as_vector([Fnum1ext, Gnum1ext]), v * n) * ds\
            + dot(as_vector([Fnum2int, Gnum2int]) - 1 / (2 * sqrt(Lambda)) * jump(C1, n), jump(CV1, n)) * dS\
            + dot(as_vector([Fnum2ext, Gnum2ext]), CV1 * n) * ds\
            + dot(as_vector([Fnum3int, Gnum3int]) - 1 / (2 * sqrt(Lambda)) * jump(C2, n), jump(CV2, n)) * dS\
            + dot(as_vector([Fnum3ext, Gnum3ext]), CV2 * n) * ds\
#putting it together
gradP = -5
F = u * v * dx - un3 * v * dx\
                 + C1 * CV1 * dx - un1 * CV1 * dx\
                 + C2 * CV2 * dx - un2 * CV2 * dx\
                 + dt * (-1 * result\
                 - (C1 - 1) * v.dx(0) * dx - 1 / Lambda * u * CV1.dx(0) * dx \
                 - (C2 - 1) * v.dx(1) * dx - 1 / Lambda * u * CV2.dx(1) * dx\
                 + gradP * v * dx + 1 / Lambda * (C1 - 1) * CV1 * dx\
                 + 1 / Lambda * (C2 - 1) * CV2 * dx)
T_end = 2
numsteps = 2000
dt = T_end / num_steps
vtkfile = File(filename)
t = 0
U1 = Function(V)
a = lhs(F)
L = rhs(F)
for n in range(num_steps):
            t += dt
            solve(a == L, U1, bc)
            print("time: ",t)
            vtkfile << (U.sub(2), t)
            u_n.assign(U)

```

I know you have to rewrite the boundary terms in terms of jumps and avgs and stuff but I suppose there lies my problem. Thank you for taking the time to read through this massive post 🙂 Every help is very much appreciated