The DistHelmholtz class¶

The DistHelmholtz<R> class is responsible for abstracting the formation and usage of the parallel sweeping preconditioner (where R is the underlying real datatype, usually double). For simplicity, this implementation uses a 7-point finite-difference stencil over a three-dimensional uniform grid. Without loss of generality (at least one face must have PML for the sweeping preconditioner to apply), PML is required on the bottom face of the domain, and the sweeps proceed from the bottom to the top of the domain. The other five boundary conditions may each be chosen to be either PML or zero Dirichlet.

type enum Boundary¶

• PML
• DIRICHLET

type struct Discretization<R>¶

R omega¶

The frequency of the Helmholtz problem in \(\frac{\mbox{rad}}{\mbox{sec}}\), where

\[\left[-\Delta - \frac{\omega^2}{c^2(x)}\right]u = f(x),\]

\(c(x)\) is the spatially varying wave propagation speed (usually in \(\frac{\mbox{km}}{\mbox{sec}}\) for seismic problems), and \(f(x)\) is the forcing function (which is implicitly equal to the physical forcing function scaled by the inverse of the spatially varying stiffness).

int nx¶

The number of grid points in the \(X\) direction.

int ny¶

The number of grid points in the \(Y\) direction.

int nz¶

The number of grid points in the \(Z\) direction.

R wx¶

The physical length of the domain in the \(X\) direction (usually in kilometers for seismic problems).

R wy¶

The physical length of the domain in the \(Y\) direction.

R wz¶

The physical length of the domain in the \(Z\) direction.

Boundary frontBC¶

Which boundary condition is to be applied on the front face of the domain (i.e., \(x=0\)).

Boundary backBC¶

Which boundary condition is to be applied on the back face of the domain (i.e., \(x=w_x\)).

Boundary leftBC¶

Which boundary condition is to be applied on the right face of the domain (i.e., \(y=0\)).

Boundary rightBC¶

Which boundary condition is to be applied on the right face of the domain (i.e., \(y=w_y\)).

Boundary topBC¶

Which boundary condition is to be applied to the top face of the domain (i.e., \(z=0\)).

int bx¶

The number of grid points of PML to be used for each boundary condition in the \(X\) direction.

int by¶

The number of grid points of PML to be used for each boundary condition in the \(Y\) direction.

int bz¶

The number of grid points of PML to be used for each boundary condition in the \(Z\) direction.

R sigmax¶

The maximum imaginary value of the complex coordinate-stretching for applying PML in the \(X\) direction. As a rule of thumb, \(1.5\, \mbox{diam}(\Omega)\) is a decent first guess.

R sigmay¶

The maximum imaginary value of the complex coordinate-stretching for applying PML in the \(Y\) direction.

R sigmaz¶

The maximum imaginary value of the complex coordinate-stretching for applying PML in the \(Y\) direction.

Discretization(R frequency, int xSize, int ySize, int zSize, R xWidth, R yWidth, R zWidth)¶

A constructor which assumes PML on all boundaries and attempts to select decent choices for the PML sizes (5 grid points) and coordinate-stretching magnitudes (\(1.5 \max\{w_x,w_y,w_z\}\)).

Discretization(R frequency, int xSize, int ySize, int zSize, R xWidth, R yWidth, R zWidth, Boundary front, Boundary right, Boundary back, Boundary left, Boundary top)¶

A constructor which attempts to select decent choices for the PML sizes (5 grid points) and coordinate-stretching magnitudes (\(1.5 \max\{w_x,w_y,w_z\}\)).

Discretization(R frequency, int xSize, int ySize, int zSize, R xWidth, R yWidth, R zWidth, Boundary front, Boundary right, Boundary back, Boundary left, Boundary top, int xPMLSize, int yPMLSize, int zPMLSize)¶

A constructor which attempts to select decent choices for the complex coordinate-stretching magnitudes (\(1.5 \max\{w_x,w_y,w_z\}\)).

Discretization(R frequency, int xSize, int ySize, int zSize, R xWidth, R yWidth, R zWidth, Boundary front, Boundary right, Boundary back, Boundary left, Boundary top, int pmlSize, double sigma)¶

A constructor which accepts a fixed value for the PML sizes and coordinate-stretching magnitudes (and applies each in all three directions).

type enum PanelScheme¶

• CLIQUE_LDL_1D: Uses triangular solves with one-dimensional frontal distributions for each of the sparse-direct subdomain solves. This is significantly slower than the next option on parallel machines.
• CLIQUE_LDL_SELINV_2D: Uses selective inversion (Raghavan et al.) for each of the subdomains in order to allow for fast subdomain solves. This is slightly less stable than CLIQUE_LDL_1D, but it increases the performance of the subdomain solves by orders of magnitude on large numbers of processes (with a small penalty in the setup time).

type class DistHelmholtz<R>¶

DistHelmholtz(const Discretization<R>& disc, mpi::Comm comm, R damping=7.5, int numPlanesPerPanel=4, int cutoff=12 )¶

Constructs a DistHelmholtz instance using the specified discretization scheme, communicator, positive imaginary shift (damping), subdomain size (numPlanesPerPanel), and maximum two-dimensional subdomain size for the analytical nested dissection (cutoff). Keep in mind that modifying damping will effect the effectiveness of the preconditioner, and numPlanesPerPanel is essentially a memory and performance tuning parameter (though it can have a minor effect on convergence rates).

void Initialize(const DistUniformGrid<R>& velocity, PanelScheme panelScheme=CLIQUE_LDL_SELINV_2D )¶

Performs the subdomain factorizations so that the preconditioner can be quickly applied.

void Solve(DistUniformGrid<Complex<R>>& B, int m=20, R relTol=1e-4, bool viewIterates=false ) const¶

Overwrites each of the right-hand sides stored in B with their approximate solution via GMRES(m) with the specified relative residual tolerance. Thus, GMRES will restart every m iterations and will continue until

\[\|A x_i - b_i\|_2 \le \mbox{relTol}\|b_i\|_2\]

for every right-hand side, \(b_i\).

void Finalize() const¶

Frees up all significant resources allocated by the class.