Solving wave equation using reduction of order

Question

I try so solve the wave equation

$$
\ddot u(x,t) - \Delta u(x,t) = f(x,t) \text{ on } D \times (0, \infty)
$$

$$
u = 0 \text{ on } \partial D \times (0, \infty)
$$
$$
u(\cdot, 0) = u_0, \qquad \dot u(\cdot, 0) = v_0 \text{ on } D
$$

use the approach

$$
\frac{\partial}{\partial t}
\begin{bmatrix} u\ v \end{bmatrix}
=
\begin{bmatrix}
0 & I \
\Delta & 0
\end{bmatrix}
\begin{bmatrix} u\ v \end{bmatrix}
+
\begin{bmatrix}
0\
f(x,t)
\end{bmatrix}
$$

I think I got most of it done already, including the weak fomrulation BUT if I'm trying to get the actual solution from my mixed vector

V_1 = FunctionSpace(mesh, 'Lagrange', 1)
V_2 = FunctionSpace(mesh, 'Lagrange', 1)
VV = MixedFunctionSpace([V_1, V_2])

u_sol = interpolate(u_ini, V_1)
v_sol = interpolate(v_ini, V_2)

for i in xrange(N_k):
    solve(A, u.vector(), b)
    tmp_u, tmp_v = u.split()
    u_sol.assign(tmp_u) # <- This is the problem
    v_sol.assign(tmp_v)

I get the error

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics@fenicsproject.org
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to initialize vector of degrees of freedom for function.
*** Reason:  Cannot re-initialize a non-empty vector. Consider creating a new function.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0
*** 
*** DOLFIN version: 1.6.0
*** Git changeset:  
*** -------------------------------------------------------------------------

What am I missing here?

Answer 1

It seems I have to use a mixed interpolation as initial condition, assign this and only aplit in the end.
Here is a complete working code.

# From http://fenicsproject.org/qa/3441/time-evolution-of-the-scalar-wave-equation
from dolfin import *
import plot

mesh = UnitSquareMesh(10, 10)
Va = FunctionSpace(mesh, "Lagrange", 1)
Vb = FunctionSpace(mesh, "Lagrange", 1)
V  = Va*Vb

t  = 0.
dt = 0.60*mesh.hmin()

def boundary(x, on_boundary):
  return on_boundary

dbc = DirichletBC(V, Constant(("0.", "0.")), boundary)

class VolumeInitialCondition(Expression):
  def eval(self, value, x):
    value[0] = exp(-(x[0]-0.5)*(x[0]-0.5) - (x[1]-0.5)*(x[1]-0.5))
    value[1] = 0

  def value_shape(self):
    return (2,)

(du, u)           = TrialFunctions(V)
(test_du, test_u) = TestFunctions(V)

w_prev = interpolate(VolumeInitialCondition(), V)
du_prev, u_prev = w_prev.split()

lhs1 = (test_u*u + test_du*du)*dx
lhs2 = dt*inner(grad(u), grad(test_du))*dx - dt*du*test_u*dx
rhs1 = test_u*u_prev*dx
rhs2 = test_du*du_prev*dx

A = assemble(lhs1 + lhs2)
b = assemble(rhs1 + rhs2)
dbc.apply(A)
dbc.apply(b)

w = Function(V)
solve(A, w.vector(), b)
output_file = File("scalar_results_euler_back.pvd", "compressed")
wave = w.split()[1]
wave.rename("WaveFunction", wave.name())
output_file << (wave, t)

for step in range(100):
  # store previous values, increment time step
    t += dt
    dut, ut = w.split(deepcopy = True)
    du_prev.assign(dut)
    u_prev.assign(ut)

  # assemble forms and solve
    A = assemble(lhs1 + lhs2)
    b = assemble(rhs1 + rhs2)
    dbc.apply(A)
    dbc.apply(b)
    solve(A, w.vector(), b)

  # save information
    wave = w.split()[1]
    wave.rename("WaveFunction", wave.name())
    output_file << (wave, t)