DG variational form

1

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:

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 :slight_smile: Every help is very much appreciated

1 Like