DOLFIN demo to solve the Helmholtz equation on a unit square

helmholtz_demo.py

"""This demo program solves Helmholtz's equation

  - div grad u(x, y) + u(x,y) = f(x, y)

on the unit square with source f given by

  f(x, y) = (1.0 + 8.0*pi**2)*cos(x[0]*2*pi)*cos(x[1]*2*pi)

and the analytical solution

  u(x, y) = cos(x[0]*2*pi)*cos(x[1]*2*pi)
"""

__author__ = "Florian Rathgeber ([email protected])"
__date__ = "2013-09-02"
__copyright__ = "Copyright (C) 2013 Florian Rathgeber"
__license__ = "GNU LGPL Version 3.0"

# Begin demo
from dolfin import *


def helmholtz(size, n):
    # Create mesh and define function space
    mesh = UnitSquare(size, size)
    V = FunctionSpace(mesh, "CG", 1)

    # Define variational problem
    lmbda = 1
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Expression("(1 + %d*pi*pi) * cos(x[0]*pi*%d) * cos(x[1]*pi*%d)" %
                   (2*n*n, n, n))
    a = (dot(grad(v), grad(u)) + lmbda*v*u)*dx
    L = f*v*dx

    # Compute solution
    A = assemble(a)
    b = assemble(L)
    x = Function(V)
    solve(A, x.vector(), b, "cg", "jacobi")

    # Analytical solution
    u = Function(V)
    u.assign(Expression("cos(x[0]*pi*%d) * cos(x[1]*pi*%d)" % (n, n)))

    # Save solution in VTK format
    #file = File("helmholtz.pvd")
    #file << x

    return x, u, x - u

if __name__ == '__main__':
    x, u, d = helmholtz(32, 2)
    # Plot solution
    plot(x, interactive=True, title="Numerical solution")
    plot(u, interactive=True, title="Analytical solution")
    plot(d, interactive=True, title="Error")

# vim:sw=4:ts=4:sts=4:et

helmholtz_mms.py

"""An MMS test for Helmholtz's equation

  - div grad u(x, y) + u(x,y) = f(x, y)

on the unit square with source f given by

  f(x, y) = (1.0 + 128.0*pi**2)*cos(x[0]*8*pi)*cos(x[1]*8*pi)

and the analytical solution

  u(x, y) = cos(x[0]*8*pi)*cos(x[1]*8*pi)

The equation is solved for different mesh resolutions and the error
norm is reported as well as the convergence order (should converge at
2nd order).
"""

__author__ = "Florian Rathgeber ([email protected])"
__date__ = "2013-09-02"
__copyright__ = "Copyright (C) 2013 Florian Rathgeber"
__license__ = "GNU LGPL Version 3.0"

# Begin demo
from dolfin import *
from numpy import asarray, log2

from helmholtz_demo import helmholtz

meshsizes = [30, 60, 120, 240]
results = [helmholtz(sz, 8) for sz in meshsizes]
errnorms = [sqrt(assemble(d*d*dx)) for x, u, d in results]

for sz, e in zip(meshsizes, errnorms):
    print sz, 'x', sz, ':', e

print "Convergence order:", log2(asarray(errnorms[:-1])/asarray(errnorms[1:]))

# vim:sw=4:ts=4:sts=4:et