Table of Contents
In this example, the scattering of an incoming plane wave between two bodies is simulated as depicted in the following figure:
 |
| Scattering of an incoming plane wave |
The mathematical formulation leads to an integral equation based on the Helmholtz kernel. Solving such an equation involves the discrete representation of the analytic operators and, if neccessary, the computation of a preconditioner. Both shall be performed using the H-matrix technique.
Problem Description
Let \(\Omega_i \in R^3, i = 1, \ldots, r\) be convex sets with smooth boundaries \(\Gamma_i = \partial \Omega_i\). Furthermore, let \(\nu = \nu_1, \nu_2, \ldots\) with \(1 \le \nu_i \le r\) be a sequence, such that \(\nu_i \ne \nu_{i+1}\) for all \(i\).
To compute: densities \(\eta_{\nu_i}\) on \(\Gamma_{\nu_i}\) defined by
\[ \eta_{\nu_i} - \mathcal{K}_{\nu_i} \eta_{\nu_i} = g_{\nu_i}, \]
with
\[ (\mathcal{K}_{\nu_i} \eta_{\nu_i})(x) := \int_{\Gamma_{\nu_i}} \frac{\partial G(x,y)}{\partial n(x)} \eta_{\nu_i}(y) dy, \quad x \in \Gamma_{\nu_i} \]
with \(G(x,y) = -\frac{1}{2\pi} \frac{e^{ik|x-y|}}{|x-y|}\), \(k > 0\) and with \(d \in R^3,\; |d| = 1,\)
\begin{eqnarray*} g_{\nu_i}(x) := \left\{ \begin{array}{ll} 2 \frac{\partial e^{ik \langle d,x \rangle }}{\partial n(x)}, & i = 1 \\ (\mathcal{K}_{\nu_{i-1}} \eta_{\nu_{i-1}})(x), \; x \in \Gamma_{\nu_i} & i > 1 \end{array} \right. \end{eqnarray*}
With a standard Galerkin discretisation using constant ansatz functions, one obtains the following matrices
\[ K^{\nu_i,\nu_i}_{ij} := \frac{1}{2\pi} \int_{\Gamma_{\nu_i}} \int_{\Gamma_{\nu_i}} e^{ik| x - y|} \frac{1-ik|x-y|}{|x-y|^3} \langle x-y, n( x ) \rangle dy dx \]
and
\[ K^{{\nu_i},{\nu_{i-1}}}_{ij} := \frac{1}{2\pi} \int_{\Gamma_{\nu_i}} \int_{\Gamma_{\nu_{i-1}}} e^{ik| x - y|} \frac{1-ik|x-y|}{|x-y|^3} \langle x-y, n( x ) \rangle dy dx \]
This leads to the following algorithm, which iterates between the surfaces \(v_i\):
- build \(K^{1,1}, K^{2,2}, K^{2,1}\) and \(K^{1,2}\);
- build initial right hand side \(b_{1,1} := g_1\);
- for \(i = 1,\cdots\)
- solve \((I - K^{1,1}) \eta_{i,1} = b_{i,1}\);
- if \(\eta_{i,1} / \eta_{i-1,1}\) converged stop
- \(b_{i,2} := K^{2,1} \eta_{i,1}\);
- solve \((I - K^{2,2}) \eta_{i,2} = b_{i,2}\);
- \(b_{i+1,1} := K^{1,2} \eta_{i,2}\);
Grid Definition
After the standard initialisation, which also imports the π§π«π¨π‘πππ types for real and complex values,
the program begins with the definition of the grid for the surfaces \(\Gamma_1\) and \(\Gamma_2\). For this, the internal grid generation of π§π«π¨π‘πππ is used to construct a spherical grid with radius 1 centered at \((0,0,0)\):
Here, "sphere-5" specifies the grid sphere with 5 refinement levels. The second command is optional. It deforms the sphere into a disc, which extends the scattering process.
Usually, a grid in π§π«π¨π‘πππ will only define triangle normals. Since this leads to constant normal directions per triangle, the approximation of the Helmholtz DLP kernel is less accurate and may even lead to false results (think of scattering between parallel planes). Therefore, the grid file should contain vertex normals, which are used to interpolate normal directions for each point on the triangle, leading to a continuous normal direction on the surface.
If vertex normals are missing, the can be trivially added under the above assumptions since each vertex coordinate equals the corresponding normal direction. You may use add_vtx_normal of TGrid for this.
Since we need two surfaces, the process is repeated, but in addition, the second disc is shifted with respect to the z direction, such that both disc will have a positive distance:
This leads to the following configuration:
 |
| Grid with two discs |
The construction of the grid ends with the definition of the corresponding function spaces. Here, we restrict ourselves to piecewise constant ansatz:
Matrix Construction
We need to construct the matrices \(K^{1,1}, K^{2,2}, K^{2,1}\) and \(K^{1,2}\). Due to definition of the grid, the matrices \(K^{1,1}\) and \(K^{2,2}\) are identical.
For the remaining matrices, cluster trees and block cluster trees are needed:
First, the matrix \(K^{1,1}\) shall be constructed. π§π«π¨π‘πππ provides the class TAcousticScatterBF, which implements the kernel for \(G(x,y)\). Not defined so far is the wavenumber \(\kappa\), which should have a value of 6. Furthermore, ACA+ is used to compute low-rank approximations with a fixed block-wise accuracy of \(10^{-4}\).
- Remarks
- Objects such as the bilinear form, the coefficient function or the low-rank approximation, which are specific for a single matrix are encapsulated in a block to prevent name clashes with the other matrices.
The actual operator in the above equation is \(I - K^{1,1}\). Here \(I\) is not the identity, but the mass matrix, as defined by the identity kernel:
Finally, \(K^{1,1}\) is updated with the mass matrix:
The construction of \(K^{1,2}\) and \(K^{2,1}\) is similar to \(K^{1,1}\). Only the function spaces and block cluster trees will change.
 |  |  |
| \(K^{1,1}\) | \(K^{1,2}\) | \(K^{2,1}\) |
Please note, that \(K^{1,1}\) is more dense than the other two matrices. This is due to the positive distance between both discs. If this is increased (or if coarsening is applied), \(K^{1,2}\) and \(K^{2,1}\) may be represented by a single low-rank matrix.
Right-Hand Side
The initial right-hand side is defined by \(g_{\nu_i}\) above, where also the the incident field \(d\) is used:
\[ g_{\nu_i}(x) = 2 \frac{\partial e^{ik \langle d,x \rangle }}{\partial n(x)} \]
Under the above setting, this is equal to:
\[ g_{\nu_i}(x) = 2 i \kappa e^{ ik \langle d,x \rangle } \langle d, n \rangle \]
and is implemented in the following RHS function:
We use TQuadBEMRHS to compute the corresponding vector of the RHS using quadrature and afterwards reorder it for H-arithmetic:
The plane wave \(d\) is chosen to be \((1,1,0)\).
Acoustic Scattering Iteration
The scattering iteration involves to solves with \(K^{1,1}\) ( \(=K^{2,2}\)). This is accelerated using a preconditioner computed by and H-LU factorisation:
We also need a vector to hold the current and the last approximation of the solution:
The variable old_ratio is needed to check convergence.
The iteration begins with solving \((I - K^{1,1}) \eta_{i,1} = b_{i,1}\), where \(b\) is the vector of the RHS:
Convergence is tested, if the ratio between the previous solution and the current solution is below some threshold. For this ratio is tested entry-wise, e.g. \(\eta_{i,1}^j / \eta_{i-1,1}^j\) for each \(j\). This is implemented in the function compare:
If this ratio does no longer change, the iteration is stopped. The threshold in this example is \(10^{-4}\):
In addition to this, the iteration shall stop if the solution \(\eta_{i,1}\) vanishes:
Otherwise, the iteration continues with the scaling of the solution vector and the computation of the RHS for \(K^{2,2}\):
Solving with \(K^{2,2}\) ( \(=K^{1,1}\)) and computing the RHS for the next iteration finished the loop and the program:
Running the simulation will produce the following output for \(\eta_{i,1}\) and \(\eta_{i,2}\). Here, the entry-wise absolute values of the complex valued solution vectors are used for the visualisation (reload page to restart animation):
 |  |  |
| \(\eta_{i,1}\) | \(\eta_{i,2}\) |
The Plain Program
- Remarks
- Within the π§π«π¨π‘πππ distribution you will find the full program with visualisation commands.
