time evolution of the scalar wave equation

Question

I'm trying to solve the scalar wave equation with homogenous Dirichlet boundary conditions and an initial state at time zero:

$$\frac{\partial^2 u}{\partial t^2} = \nabla^2 u,$$

by turning it into two, coupled, first-order equations - and solve them using the Crank-Nicolson method:

$$\frac{u^k - u^{k -1}}{\delta t} = (w^k + w^{k -1})/2$$
$$\frac{w^k - w^{k - 1}}{\delta t} = (\nabla^2 u^k + \nabla^2 u^{k-1})/2$$

The code below does show that the volume initial condition is applied, but I've done something that causes the solution to never change … :(

from dolfin import *

CRANK_NICOLSON       = 0.5
theta = Constant(CRANK_NICOLSON)

mesh = UnitCubeMesh(10, 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):
    r = sqrt((x[0] - 0.5)**2 + (x[1] - 0.5)**2 + (x[2] - 0.5)**2)
    value[0] = exp(-(r/0.2)**2)
    value[1] = -2/0.2*r*exp(-(r/0.2)**2)

  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()

Dt = Constant(dt/2.)
lhs1 = test_u*(u - Dt*du)*dx
rhs1 = test_u*(u_prev + Dt*du_prev)*dx
lhs2 = test_du*du*dx + Dt*inner(grad(u), grad(test_du))*dx
rhs2 = test_du*du_prev*dx - Dt*inner(grad(test_du), grad(u_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.pvd", "compressed")
wave = w.split()[1]
wave.rename("WaveFunction", wave.name())
output_file << (wave, t)

for step in range(20):
  # store previous values, increment time step
  t += dt
  w_prev.assign(w)
  du_prev, u_prev = w_prev.split()

  # 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)

Comments

Hi,

I'm a bit confused by your code, since I'm used to handling wave equations a bit differently:

1) When splitting the second order wave equation up into two first order equations, I usually do this as follows:
$$
\dot{u} = \nabla \cdot v
$$
$$
\dot{v} = \nabla v
$$
That way the system is also first order with respect to the spatial derivatives.
2) You're working on a 3D mesh, but your unknown fields are scalars. This usually means you assume two of the three components of the unknown vector fields are zero, so you can compute only the non-zero component. In that case also the grad should be reduced to a single derivative, e.g., .dx(0) for the x-direction
3) Don't apply Dirichlet boundary conditions to both unknowns.

That being said, I'm not sure what's actually wrong with your code yet...
I'll try to write my own code for this simple example...

---

Thanks for the reply - my laptop is an ancient Windows machine, so I'll have to wait till Monday to try your idea out. :(

I really wanted to solve a vector wave equation, but figured that the scalar wave equation would be "easier" to start out with. Is there a reference that describes your divergence trick?

---

Hello stevenvdk and timm, I am also trying to model the seismic wave equation on a 2D mesh. You seem to have excellent previous experience with wave modelling. Can you guys pls. help me with my question below? I' ll be really grateful for your help and guidance!

http://fenicsproject.org/qa/8352/best-method-for-discretization-of-wave-equation

Answer 1

Hi,

you can try to replace

w_prev.assign(w)
du_prev, u_prev = w_prev.split()

by

dut, ut = w.split(deepcopy = True)
du_prev.assign(dut)
u_prev.assign(ut)

That will make your results "propagate". However, since the rest of your code needs some work (see my comments), the result is not yet correct...

Regards

Comments

Thanks - sometimes the Python script makes everything a bit mysterious because I probably really need to understand what is going on with the data structures in the background. Would I be better off coding this in C++?

---

Hey! that got it to propagate the results … thanks for the advice on the change.