I am having trouble specifying a Dirichlet boundary condition when I attempt to solve the following problem using a discontinuous Galerkin approach adapted from (link):
\frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x} - a\frac{\partial^2u}{\partial x^2} = 0
on a 1D interval mesh from (0,2\pi) with a and c = 1.
I know that under periodic boundary conditions with an initial condition of \mathrm{sin}(x), the solution should be u(x,t) = e^{-at}\mathrm{sin}(x-ct). As FEniCS currently doesn’t support periodic boundary conditions for DG elements, I am attempting to use the exact solution as a Dirichlet boundary condition to check if the DG approximation is correct. However, it doesn’t appear that I am able to define a Dirchlet boundary condition that alters the solution (i.e. even setting the Dirichlet boundary to be a constant doesn’t alter the boundary values of the solution).
I am unsure if this is due to the way I have defined my variational problem, or perhaps I need to use a different approach to properly implement the Dirichlet condition.
Any insights you can offer would be greatly appreciated.
I am using DOLFIN 2019.1.0 on Ubuntu 18.04
%matplotlib inline
import time
import os
import math
from dolfin import *
import matplotlib.pyplot as plt
# u_t + c*u_x - a*u_xx = 0
# Parameters
parameters['form_compiler']['cpp_optimize'] = True
parameters['form_compiler']['optimize'] = True
parameters["ghost_mode"] = "shared_facet"
t = 0
t_end = 2
dt = 0.05
# Create mesh
mesh = IntervalMesh(40, 0, 2*pi)
# Define function spaces
V = FunctionSpace(mesh, "DG", 1)
# Set initial condition and define c
ic= Expression('sin(x[0])', degree=6, domain=mesh)
c = Expression('1', degree = 0, domain = mesh)
# Define unknown and test function(s)
v = TestFunction(V)
u = TrialFunction(V)
# Previous Solution
u0 = Function(V)
u0 = interpolate(ic,V )
# Penalty
h = CellDiameter(mesh)
n = FacetNormal(mesh)
alpha = Constant(1)
# Upwinding
cv = as_vector((c,))
cn = (dot(cv, n) + abs(dot(cv, n)))/2.0
# Boundary conditions (using exact solution instead of periodic)
b_value = Expression('exp(-t)*sin(x[0]-t)', t=t, degree=2)
bc = DirichletBC(V, b_value, 'on_boundary')
F = (1/dt)*(u-u0)*v*dx + dot(grad(v), grad(u) -cv*u)*dx\
+ (alpha/avg(h))*dot(jump(u, n), jump(v, n))*dS \
- dot(avg(grad(u)), jump(v, n))*dS \
- dot(jump(u, n), avg(grad(v)))*dS \
+ (jump(v))*(cn('+')*u('+') - cn('-')*u('-'))*dS + v*cn*u*ds
u = Function(V)
sys_l, sys_r = lhs(F), rhs(F)
solve(sys_l==sys_r, u, bc)
u.assign(u0)
u.rename("u", "u")
# Time-stepping
t = 0.0
while t < t_end:
print("t =", t, "end t=", t_end)
# Compute
b_value.t = t
bc = DirichletBC(V, b_value, 'on_boundary')
sys_l, sys_r = lhs(F), rhs(F)
solve(sys_l==sys_r, u, bc)
plot(u, label='DG')
u_e = Expression('exp(-t)*sin(x[0]-t)', degree=6, domain=mesh, t=t)
u_sol = project(u_e, V)
plot(u_sol, label='Exact')
plt.legend()
plt.ylim(-1,1)
plt.show()
# Move to next time step
u0.assign(u)
t += dt
