HLIBpro 3.1
|
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.
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}{4\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} 4 \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}{4\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}{4\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\):
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.
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:
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 THelmholtzDLPBF
, which implements the standard and adjoint Helmholtz DLP kernel, for which the latter is used for \(G(x,y)\). Not defined so far is the wavenumber \(\kappa\), which should have a value of 16. Furthermore, ACA+ is used to compute low-rank approximations with a fixed block-wise accuracy of \(10^{-4}\).
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.
The initial right-hand side is defined by \(g_{\nu_i}\) above, where also the incident field \(d\) is used:
\[ g_{\nu_i}(x) = 4 \frac{\partial e^{ik \langle d,x \rangle }}{\partial n(x)} \]
Under the above setting, this is equal to:
\[ g_{\nu_i}(x) = 4 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)\).
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}\) |