Helmholtz Problem
The partial differential equation to be solved for
has the form
 |
(21) |
and we set
and  |
(22) |
With
the radiation condition defined by Equation (1.21)
takes the form
on  |
(23) |
The partial differential Equation (1.22) together with boundary conditions of Equation (1.24)
is called the Helmholtz equation .
We want to use the LinearPDE class provided by esys.escript to define and solve a general linear,steady, second order PDE such as the
Helmholtz equation. For a single PDE the LinearPDE class supports the following form:
 |
(24) |
where we show only the coefficients relevant for the problem discussed here. For the general form of
single PDE see Equation (4.1).
The coefficients
, and
have to be specified through Data objects in the
general FunctionSpace on the PDE or objects that can be converted into such Data objects.
is a rank-2 Data object and
and
are scalar.
The following natural
boundary conditions are considered on
:
 |
(25) |
Notice that the coefficient
is the same like in the PDE Equation (1.25).
The coefficients
and
are
each a scalar Data object in the boundary FunctionSpace. Constraints for the solution prescribing the value of the
solution at certain locations in the domain. They have the form
where  |
(26) |
and
are each scalar Data object where
is the characteristic function
defining where the constraint is applied.
The constraints defined by Equation (
) override any other condition set by
Equation (1.25) or Equation (1.26).
The Poisson class of the esys.escript.linearPDEs module,
which we have already used in Chapter 1.2, is in fact a subclass of the more general
LinearPDE class. The esys.escript.linearPDEs module provides a Helmholtz class but
we will make direct use of the general LinearPDE class.
By inspecting the Helmholtz equation
(1.22) and boundary condition (1.24) and
substituting
for
we can easily assign values to the coefficients in the
general PDE of the LinearPDE class:
 |
(27) |
is the Kronecker symbol defined by
for
and 0 otherwise. Undefined coefficients are assumed to be not present.
In this diffusion example we do not need to define a characteristic function
because the
boundary conditions we consider in Equation (1.24) are just the natural boundary
conditions which are already defined in the LinearPDE class (shown in Equation (1.26)).
The Helmholtz equation can be set up by following way
:
where we assume that mydomain is a Domain object and
kappa omega eta and g are given scalar values
typically float or Data objects. The setValue method
assigns values to the coefficients of the general PDE. The getSolution method solves
the PDE and returns the solution u of the PDE. kronecker is esys.escript function
returning the Kronecker symbol.
The coefficients can set by several calls of setValue where the order can be chosen arbitrarily.
If a value is assigned to a coefficient several times, the last assigned value is used when
the solution is calculated:
In some cases the solver of the PDE can make use of the fact that the PDE is symmetric where the
PDE is called symmetric if
 |
(28) |
Note that
and
may have any value and the right hand sides
,
as well as the constraints
are not relevant. The Helmholtz problem is symmetric.
The LinearPDE class provides the method checkSymmetry to check if the given PDE is symmetric.
Unfortunately, a checkSymmetry is very expensive and is recommendable to use for
testing and debugging purposes only. The setSymmetryOn method is used to
declare a PDE symmetric:
Now we want to see how we actually solve the Helmholtz equation.
on a rectangular domain
of length
and height
. We choose a simple test solution such that we
can verify the returned solution against the exact answer. Actually, we
take
(here
) and then calculate the right hand side terms
and
such that
the test solution becomes the solution of the problem. If we assume
as being constant,
an easy calculation shows that we have to choose
. On the boundary we get
.
So we have to set
. The following script helmholtz.py
which is available in the example directory
implements this test problem using the esys.finley PDE solver:
To visualize the solution `x0. xml' just use the command
and it is easy to see that the solution
is calculated.
The script is similar to the script poisson.py discussed in Chapter 1.2.
mydomain.getNormal() returns the outer normal field on the surface of the domain. The function Lsup
imported by the from esys.escript import * statement and returns the maximum absolute value of its argument.
The error shown by the print statement should be in the order of
. As piecewise bi-linear interpolation is
used by esys.finley approximate the solution and our solution is a linear function of the spatial coordinates one might
expect that the error would be zero or in the order of machine precision (typically
).
However most PDE packages use an iterative solver which is terminated
when a given tolerance has been reached. The default tolerance is
. This value can be altered by using the
setTolerance of the LinearPDE class.
esys@esscc.uq.edu.au