Helmholtz Problem

The partial differential equation to be solved for $ T^{(n)}$ has the form

$\displaystyle \omega T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f$ (21)

and we set

$\displaystyle \omega=\frac{\rho c\hackscore p}{h}$    and $\displaystyle f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.$ (22)

With $ g=\eta T\hackscore{ref}$ the radiation condition defined by Equation (1.21) takes the form

$\displaystyle \kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g - \eta T^{(n)}$ on $\displaystyle \Gamma$ (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:

$\displaystyle -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u = Y \; .$ (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 $ A$, and $ Y$ have to be specified through Data objects in the general FunctionSpace on the PDE or objects that can be converted into such Data objects. $ A$ is a rank-2 Data object and $ D$ and $ Y$ are scalar. The following natural boundary conditions are considered on $ \Gamma$:

$\displaystyle n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+d u= y \;.$ (25)

Notice that the coefficient $ A$ is the same like in the PDE Equation (1.25). The coefficients $ d$ and $ y$ 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

$\displaystyle u=r$    where $\displaystyle q>0$ (26)

$ r$ and $ q$ are each scalar Data object where $ q$ 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 $ u$ for $ T^{(n)}$ we can easily assign values to the coefficients in the general PDE of the LinearPDE class:

\begin{displaymath}\begin{array}{llllll} A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\ d=\eta & y= g & \\ \end{array}\end{displaymath} (27)

$ \delta\hackscore{ij}$ is the Kronecker symbol defined by $ \delta\hackscore{ij}=1$ for $ i=j$ and 0 otherwise. Undefined coefficients are assumed to be not present.[*] In this diffusion example we do not need to define a characteristic function $ q$ 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 [*] :
\begin{python}
mypde=LinearPDE(mydomain)
mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g)
u=mypde.getSolution()
\end{python}
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:
\begin{python}
mypde=LinearPDE(mydomain)
mypde.setValue(A=kappa*kronecker(mydoma...
...g)
mypde.setValue(d=2*eta)  ...
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

$\displaystyle A\hackscore{jl}=A\hackscore{lj}\;.$ (28)

Note that $ D$ and $ d$ may have any value and the right hand sides $ Y$, $ y$ 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.
\begin{python}
mypde=LinearPDE(mydomain)
mypde.setValue(A=kappa*kronecker(mydoma...
...=kronecker(mydomain)[0])
print mypde.checkSymmetry()  ...
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:
\begin{python}
mypde = LinearPDE(mydomain)
mypde.setValue(A=kappa*kronecker(mydomain))
mypde.setSymmetryOn()
\end{python}
Now we want to see how we actually solve the Helmholtz equation. on a rectangular domain of length $ l\hackscore{0}=5$ and height $ l\hackscore{1}=1$. We choose a simple test solution such that we can verify the returned solution against the exact answer. Actually, we take $ T=x\hackscore{0}$ (here $ q\hackscore H = 0$) and then calculate the right hand side terms $ f$ and $ g$ such that the test solution becomes the solution of the problem. If we assume $ \kappa$ as being constant, an easy calculation shows that we have to choose $ f=\omega \cdot x\hackscore{0}$. On the boundary we get $ \kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. So we have to set $ g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script helmholtz.py which is available in the example directory implements this test problem using the esys.finley PDE solver:
\begin{python}
from esys.escript import *
from esys.escript.linearPDEs import Li...
...lution()
print ''error is '',Lsup(u-x[0])
saveVTK(''x0.xml'',sol=u)
\end{python}
To visualize the solution `x0. xml' just use the command
\begin{python}
mayavi -d u.xml -m SurfaceMap &
\end{python}
and it is easy to see that the solution $ T=x\hackscore{0}$ 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 $ 10^{-7}$. 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 $ \approx 10^{-15}$). However most PDE packages use an iterative solver which is terminated when a given tolerance has been reached. The default tolerance is $ 10^{-8}$. This value can be altered by using the setTolerance of the LinearPDE class.

esys@esscc.uq.edu.au