Direct structured solvers and time-harmonic waves
In the invere boundary value problem, to probe the Dirichlet-to-Neumann map, one needs to solve the system of equations describing time-harmonic waves for a large number of right-hand sides. We use (compact) FD (with mass lumping) and CG-FEM methods for the discretization. A typical region comprises a three-dimensional structured or unstructured tetrahedral mesh consisting of about 109 elements.
We introduced massively parallel structured direct solvers for this purpose using Hierarchically Semi-Separable matrices, transforming the way time-harmonic waves are computed approaching linear complexity. We exploit the fact that data are acquired with limited accuracy and use techniques to control finite accuracy computation. Our methods involve multiple layers of hierarchical tree structures and advanced graph and sparse matrix techniques.
- Hierarchically Semi-Separable (HSS) matrices, structured multifrontal solvers. HSS matrix algorithms have been emerging in the construction of the superfast direct solvers for both dense and sparse linear systems. We developed a set of new parallel algorithms for the key HSS operations that are used for solving such large systems. These include the parallel rank-revealing QR factorization, the HSS constructions with hierarchical compression, the ULV HSS factorization, and the HSS solutions [105]. We coupled these to the multifrontal method, using a matrix re-ordering following nested dissection based domain decomposition, and designed and implemented massively parallel algorithms [89,94,100]. The rank patterns of the off-diagonal blocks are precisely such that we obtain a storage requirement which is nearly O(N) (N signifies the size of the matrix) and a computational complexity of the factorization which falls just below O(N 4/3) in dimension three.
- The method in [105] follows a geometric approach with fits discretizations based on rectangular grids or structured meshes. We have now completed the development of an algebraic approach, which is amenable to the use of unstructured tetrahedral meshes (arising in the conditional Lipschitz stability estimates for the inverse boundary value problem and multi-level schemes).
The direct structured solvers also have applications in certain problems of computational geodynamics; the un- structured tetrahedral meshes – with multi-scale refinement and deformation – provide an interface between these and (seismic) inverse problems.
We developed distributed memory parallel randomized HSS algorithms, which avoid the storage of any dense matrix. The parallelism in the HSS tree structure is exploited by overlapping the computations on each separate subtree using different subgroups of processors. Within each subgroup, a 2D process grid is created for the block-cyclic matrix distribution of each HSS generator. If r is the rank bound of the off-diagonal blocks and n is the matrix size, the parallel run time for randomized HSS algorithms is essentially O(r) with O(rn) processors.
- Nonlinear eigenvalue problems and normal modes. We anticipate that the structured matrix approach will lead to a superfast algorithm to solve nonlinear eigenvalue problems (three-dimensional anisotropic, rotating, self- gravitating Earth models) in prescribed frequency intervals. This would make it feasible to compute directly normal modes in a three-dimensional heterogeneous earth.
The structured multifrontal factorization techniques apply to numerically solving certain integral equations, such as the Lippman-Schwinger equation and, via a Neumann series, the integral equation driving (time-harmonic) surface-related multiple elimination.
Time-domain solvers
In the time domain, we developed a novel Discontinuous Galerkin (DG) method for acousto-elastic waves allowing general anisotropy. In my program, the use of the DG method is primarily motivated by the introduction of unstructured tetrahedral meshes (and polyhedral interfaces) through the relevant hyperbolic inverse problems subject to conditions of well-posedness. This method yields a block-diagonal mass matrix.
We started from a coupled first-order elastic strain-velocity and acoustic velocity-pressure formulation. We introduced a novel penalty flux and avoided diagonalization. Our penalty flux accommodates fluid-solid boundaries. We obtained a convergence result. In the implementation, we made use of a fast implicit-explicit time integration approach so as to handle a local multi-scale refinement of the mesh, again, arising in iterative schemes for hyperbolic inverse problems. (We note that the structured multifrontal factorization using HSS matrices directly applies to the fast computation of the propagator, that is, a complex exponential map, associated the half scalar-wave equation.)
We have designed a parallel computational platform from the laboratory (microstructures such as fracture systems) and exploration (geological structures, salt tectonics) to global (planetary) scale.
Model structures, flexible updating and inverse problems
It seems to be natural to integrate (the modelling of) structural geological and geodynamical processes with domain partitioning (meshing) and associated sparse multi-scale model representations for which the mentioned Lipschitz stability estimates hold. To constrain the deformation during updating in iterative methods, we exploit connections with deforming meshes (and mesh refinement). Techniques based on topological completion formulas obtained by a homology criterion are applicable.
