Hello all. I am new to both fenics and finite element problems in general. I am trying to replicate the results found here: Bitbucket using implicit Euler time-stepping rather than spacetime finite elements.
However, my code as constructed causes the appearance of an artefact near the shock that is resistant to changes in the size of the timestep.
I’m unsure if this problem is caused by the change in approach to the time dimension or an incorrect specification of the Lax-Friedrichs flux
import matplotlib.pyplot as plt
import ufl
from dolfin import (
Constant, IntervalMesh, SubDomain, DOLFIN_EPS, project, inner,
as_vector, PETScSNESSolver, dot, dS, NonlinearProblem, plot,
FunctionSpace, TrialFunction, TestFunction, Expression, MeshFunction,
Measure, FacetNormal, derivative, assemble, dx, interpolate, Function,
solve, DirichletBC, avg, jump, nabla_grad, grad)
import time
import numpy as np
start = time.time()
#domain parameters
x_R = 1.
x_L = 0.
N_el = 64
L = x_R - x_L
cell_size = (x_R-x_L)/N_el
#Parameters for time marching algorithm
T = 0.5
dt = 0.5/64
dt_ufl = Constant(dt)
#Viscosity
nu = Constant(0.01)
# Mesh and function space.
mesh = IntervalMesh(N_el, x_L, x_R) # mesh
V = FunctionSpace(mesh, 'DG', 1) #function space
# Boundary and initial conditions
# Set up IC
u0_expr = Expression('max(sin(2.0*pi*x[0] ), 0.0 )', element=V.ufl_element())
Dirichelet_inlet = Constant(0) # Dirichlet BC at inlet
#functions to define the boundary subdomains
class Inlet_BC(SubDomain):
def inside(self, x, on):
return (abs(x[0] - x_L) < DOLFIN_EPS) and on
class Outlet_BC(SubDomain):
def inside(self, x, on):
return (abs(x[0] - x_R) < DOLFIN_EPS) and on
class Burgers_Variational_form():
def __init__(self, x_L, x_R, N_el, dt, nu, mesh, V, u0, Inlet_BC):
self.x_L = x_L
self.x_R = x_R
self.N_el = N_el
self.dt = dt
self.nu = nu
self.dt_ufl = Constant(dt)
# Mesh and function space.
self.mesh = mesh
self.V = V
self.u = Function(V)
self.phi = TestFunction(V)
self.n = FacetNormal(mesh)
#BCs and ICs
self.u_old = interpolate(u0, V)
# self.u.assign(self.u_old)
# Boundary conditions
exterior_bdries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
inlet_ds = Inlet_BC()
inlet_ds.mark(exterior_bdries, 1)
outlet_ds = Outlet_BC()
outlet_ds.mark(exterior_bdries, 2)
self.ds = Measure('ds', domain=mesh, subdomain_data=exterior_bdries)
def Fc(self, u):
flux = u**2/2
return flux
# Define local-Lax Friedrichs flux for the convective term
def numerical_cflux(self, u_p, u_m, n):
# note that n^+ = -n^-
C = ufl.Max(abs(u_p), abs(u_m))
# flux = {{fc(u)}} + C/2*[[u]]
flux = (self.Fc(u_p) + self.Fc(u_m))/2 + (C/2)*(u_p*n[0] - u_m*n[0])
return flux
def domain_numerical_flux(self, u, n):
interior_numerical_fluxes = inner(self.numerical_cflux(u('+'), u('-'), n('+')), self.phi('+') - self.phi('-'))*dS
# Exterior Dirichlet condition
inlet_flux_BC = self.numerical_cflux(u, Dirichelet_inlet, n)*self.phi*self.ds(1)
# Exterior Neumann condition at outlet
outlet_flux_BC = self.numerical_cflux(u, u, n)*self.phi*self.ds(2)
return interior_numerical_fluxes + inlet_flux_BC + outlet_flux_BC
#volumetric term
def volume_term(self, u, u_old):
volume = (inner(u - u_old, self.phi) - self.dt_ufl * inner(self.Fc(u), self.phi.dx(0) )) * dx
return volume
def get_residual(self):
residual = self.volume_term(self.u, self.u_old) + self.dt_ufl * self.domain_numerical_flux(self.u, self.n)
J = derivative(residual, self.u)
return residual, J
if __name__ == '__main__':
form = Burgers_Variational_form(x_L, x_R, N_el, dt, nu, mesh, V, u0_expr, Inlet_BC)
residual, J = form.get_residual()
u = form.u
u_old = form.u_old
class BurgersProblem(NonlinearProblem):
def F(self, b, x):
assemble(residual, tensor=b)
def J(self, A, x):
assemble(J, tensor=A)
burgers = BurgersProblem()
solver = PETScSNESSolver('newtonls')
# Time stepping
t = 0
i = 0
while t < T - 0.5*dt:
solver.solve(BurgersProblem(), u.vector())
t += dt
i += 1
u_old.assign(u)
plot(u)
plt.show()
A comparison between the example solution and my solution can be seen below:
Any thoughts or advice would be greatly appreciated. Thank you!