Need help to understand this working code : wave equation + crank nicolson

Question

Hello everyone,

I found this code on wave equation using Crank Nicolson and am trying to understand how the CN scheme is being implemented here. So, the form of crank nicolson scheme, that I am aware of is: where theta = 1/2 and 'k' is the time step.

$$u_t = v, v_t = \Delta u + f$$

Assembling the matrices and writing in term of Mass matrix 'M' and stiffness matrix 'A'.

The core of the code is shown below:

m = u*v*dx
M = assemble(m)
a = inner(grad(u),grad(v))*dx
A = assemble(a)
B = M + k**2*A/4.            
C = 2.*M - k**2*A/2.
D = -1.*M - k**2*A/4.

# solve linear variational problem
def solveprob(m1, m2, u, bc):
  prob = LinearVariationalProblem(m1, m2, u, bcs=bc)
  solver = LinearVariationalSolver(prob)
  solver.solve()

# Modify B to enforce boundary conditions
bc = DirichletBC(V, g, boundary)
bc.apply(B)

# compute LU factorization of B
Bsolve = LUSolver(B)
Bsolve.parameters["reuse_factorization"] = True

# initial data is L^2-projection
g.t = t0 # evaluate g at time t0
uold = Function(V)
solveprob(m, u0*v*dx, uold, bc)
vold = Function(V)
solveprob(m, v0*v*dx, vold, bc)

f.t = t0
foldvec = assemble(f*v*dx)
if saveframes: vtkfile << uold

# approximate u(k) by u0 + k v0 + k^2/2 (Lap u + f)
t = t0 + k
ucur = Function(V)
g.t = t # evaluate g at time t0 + k
F = (uold + k*vold + k**2/2.*f)*v*dx - k**2/2.*inner(grad(uold),grad(v))*dx
solveprob(m, F, ucur, bc)
f.t = t
fcurvec = assemble(f*v*dx)
if saveframes: vtkfile << ucur

# Compute solution
unew = Function(V)
tstep = 1
print 'Beginning timestepping... '
while t + k <= T:
    # Crank-Nicolson timestepping applied to u_t = v, v_t = Lap u + f, and then v eliminated:
    # (u_{n+1} - 2 u_n + u_{n-1})/k^2 = (Lap u_{n+1} + 2 Lap u_n + Lap u_{n-1})/4 + (f_{n+1} + 2 f_n + f_{n-1})/4
    t += k
    tstep += 1
    f.t = t
    fnewvec = assemble(f*v*dx)
    L = C*ucur.vector() + D*uold.vector() + k**2 * (foldvec + 2.*fcurvec + fnewvec)/4.
    g.t = t # evaluate g at new time
    bc.apply(L)
    Bsolve.solve(unew.vector(), L)
    if saveframes: vtkfile << unew
    uold.assign(ucur)
    ucur.assign(unew)
    foldvec = fcurvec
    fcurvec = fnewvec
    if tstep%10 == 0:
      print "timestep: %g" % tstep

if not saveframes: vtkfile << ucur
print 'L2 norm at final time %.4e' % sqrt(abs(assemble(unew**2*dx)))

Now, I have the following questions:

1. It is mentioned in the code that the form being implemented is:

Crank-Nicolson timestepping applied to u_t = v, v_t = Lap u + f, and then v eliminated:
$$(u_{n+1} - 2 u_n + u_{n-1})/k^2 = (\Delta u_{n+1} + 2 \Delta u_n + \Delta u_{n-1})/4 + (f_{n+1} + 2 f_n + f_{n-1})/4$$

How is the above form obtained?

2. What is the deal here with saveframes command used with 'if' statement for uold, ucur and unew? I have no idea what that is doing.

3. Where does the following come from ?

 F = (uold + k*vold + k**2/2.*f)*v*dx - k**2/2.*inner(grad(uold),grad(v))*dx
 L = C*ucur.vector() + D*uold.vector() + k**2 * (foldvec + 2.*fcurvec + fnewvec)/4.