Setting up variational formulation for DG 1 in time

Question

I'm trying to set up the variational formulation for DG1 in time for the linear test equation
$$ \frac{d}{dt}u(t) = -u(t), $$
which would be
$$\int_{I_n} \frac{d}{dt}u^n(t)v(t)dt + (u^n(t_n) - u^{n-1}(t_n))v(t_n) + \frac{\Delta t_n}{2}(u^n(t_n) v(t_n) + u^n(t_{n+1})v(t_{n+1})),$$
with $$ I_n := (t_n, t_{n+1}) $$ being the interval for the current timestep, $$u^n(t)$$ the solution associated to $$I_n$$ and $$u^{n-1}(t)$$ corresponding to the previous timestep/initial value.

I figure I can implement the derivative term by

inner(v, u_new.dx(0))*dx

and the term $$\frac{\Delta t_n}{2}(u^n(t_n) v(t_n) + u^n(t_{n+1})v(t_{n+1}))$$ by

0.5*timestep*inner(v, u_new)*ds

but I'm a little clueless on how to implement the remaining term. I don't quite understand how to use the restrictions like v('+'), v('-') to make this work.

Minimal example until now:

mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh, 'DG', 1)
u_old = Function(V)
u_old.interpolate(Expression('1'))
v = TestFunction(V)
u_new = Function(V)
timestep = Constant(1)
F = inner(v, u_new.dx(0))*dx + 0.5*timestep*inner(v, u_new)*ds 
solve(F == 0, u_new)
print('result', u_new.vector().array()) # (0.4, 0.8) expected

Any help/advice would be welcome!

Comments

Could solve the problem by using Subdomains to integrate over the specific facets.

---

Hello Peter, would you mind posting your solution? I'm intested in a similar case... thanks!

Answer 1

Solution I found:

T = 1
n = 100
timestep = Constant(T/n)

mesh = UnitIntervalMesh(1)
V = FunctionSpace(mesh, 'DG', 1)
class bound_left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)
class bound_right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

leftboundary = bound_left()
rightboundary = bound_right()

boundaries = FacetFunction('size_t', mesh)
leftboundary.mark(boundaries, 1)
rightboundary.mark(boundaries, 2)
ds = Measure('ds')[boundaries]

u_new, u_old = Function(V), Function(V)
u_old.interpolate(Expression('1'))
v = TestFunction(V)
F = (inner(u_new.dx(0), v)*dx + 0.5*timestep*inner(u_new, v)*ds + (u_new - u_old)*v*ds(1))
for i in range(n):
    solve(F == 0, u_new)
    u_old.interpolate(Constant(assemble(u_new*ds(2))))

I use subdomains to be able to evaluate a given Function at a specific side of my intervals via ds(1), ds(2). I resolved the problem of evaluating $$ u^{n-1}(t_n)v(t_n)$$, which needs right side for u and left side for v, by interpolating u_old to a constant function with the right value after every step.