Unstructured Meshes, Deformation, Multiscale Refinement

Fast algorithms, massively parallel computing

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.

Discontinuous Galerkin Method: SEAM 3D

Fast algorithms, massively parallel computing

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.

Massively Parallel Structured Direct Helmholtz Solver

Fast algorithms, massively parallel computing

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.

Conditional Lipschitz Stability, Sparsity, Convergence Radius

Nonlinear inverse (boundary value) problems: Conditional well-posedness, iterative methods

Part of my research program has been devoted to developing new, comprehensive quantitative probes of our planet’s interior which include novel convergent nonlinear iterative techniques (surpassing so-called standard ‘adjoint state’ methods). One of our most recent breakthroughs concerns the understanding of the fundamental role of conditional Hölder (or Lipschitz) stability estimates of the nonlinear seismic inverse problem.

Novel iterative methods, conditional Lipschitz stability and convergence

Forward modelling and optimization to match waveforms is an approach that was introduced in seismology more than four decades ago. Even though this approach, commonly using implementations based on the ‘adjoint state’ method, is being rather widely applied today, no results on the recovery of coefficients of models, however, were obtained. The first fundamental question to be addressed is whether the inverse boundary value problem corresponding with idealized seismic data has a unique solution. The second question is whether explicit reconstructions can be obtained, and the third question is what one can say about the case of partial data, also subject to sampling. The fourth question concerns the incorporation of inaccurate data and uncertainty quantification. All these questions are intimately interconnected. The second question typically involves iterative methods (‘full waveform inversion’) and leads to a local result.

We succeeded in obtaining conditional convergence of iterative methods results; we obtained explicit expressions for the radii of convergence [96]. The convergence relies on a Hölder, or Lipschitz, stability estimate for (essentially well-posedness of) the inverse problem. These insights have a significant impact on our understanding and develop- ment of optimization techniques very different from optimization with traditional regularization strategies.

  • Lipschitz stability estimates: scalar waves. The time-harmonic or ‘frequency-domain’ formulation lends itself to obtain conditional Lipschitz stability estimates for the nonlinear seismic inverse boundary value problem, even assuming partial boundary data [103]. The data are the Dirichlet-to-Neumann map. Uniqueness of recovery for this inverse problem in the case of scalar waves has been known since the late 1980s under minimal regularity assumptions. We established these estimates in the case of piecewise constant models with underlying domain partitionings containing discontinuities.

    The Lipschitz stability constant grows exponentially with the number of subdomains in the domain partitioning. The stability constant appears in the expressions for the radii of convergence, in a way such that the radius tends to zero as this constant becomes large. We obtained a sharp estimate of the stability constant in terms of the number of subdomains in the domain partitioning [preprint 7]. This estimate plays a critical role in the development of a multi-level (multi-scale) scheme.

    Moreover, we obtained the Lipschitz stable recovery of polyhedral interfaces and – if the domain partitioning in the above is an (unstructured) tetrahedral mesh – the Lipschitz stable recovery of the number of elements in and the shape of the mesh [preprint 8]. Thus we fundamentally relaxed knowledge about the parametrization.

    The Dirichlet-to-Neumann map is related to a single layer potential operator, the kernel of which is directly related to seismic data.

  • Time-harmonic elastic waves. We obtained the unique recovery of Lame´ parameters and density that are piece- wise constant on an unstructured tetrahedral mesh from the Neumann-to-Dirichlet map (which corresponds with vibroseis data). We also obtained Lipschitz stability.
  • Multi-level, multi-scale schemes. Based on the conditional stability estimates we introduce convex subsets, fol- lowing a refinement of domain partitionings. We designed an iterative method which converges to the point in each subset closest to the unique solution. The distance, or approximation error, depends on frequency, and an estimate can be obtained. Mitigating the growth of the stability constants with the approximation errors controlled by frequency, we obtained a multi-level (nonlinear projected steepest descent iteration) scheme with precise conditions for convergence [109] defining a nonlinear hierarchical compressed reconstruction. Through the compression and corresponding approximation we incorporate inaccuracies in the data. Other more ad- vanced methods, in particular the Gauss-Newton method, can also be adapted for this purpose.
  • Attenuation. The Lipschitz stability estimates hold for complex (frequency dependent) coefficients and, hence, the time-harmonic formulation allows in principle to obtain localized information about the (reciprocal) Q factor, and spatial variability in physical dispersion. This might open the way to the further studying of attenuation in localized regions in the mantle, such as in the asthenosphere.

Hyperbolic inverse boundary value problem

The time-domain formulation of the inverse boundary value problem for scalar and elastic waves involves different techniques and strategies. Naturally, in this formulation, one can introduce time windows identifying structure in and emphasizing parts of the data. We have been studying this problem for coefficients with structure. On the one hand, in the case of piecewise smooth coefficients with interfaces of unknown shapes, we use connections with geometric inverse problems with corresponding conditions; on the other hand, we establish explicit connections with the time- harmonic inverse boundary value problem, a multi-level strategy, and conditional high-frequency stability estimates.

Conditional Lipschitz Stability, Sparsity, Convergence Radius

Nonlinear inverse (boundary value) problems: Conditional well-posedness, iterative methods

Part of my research program has been devoted to developing new, comprehensive quantitative probes of our planet’s interior which include novel convergent nonlinear iterative techniques (surpassing so-called standard ‘adjoint state’ methods). One of our most recent breakthroughs concerns the understanding of the fundamental role of conditional Hölder (or Lipschitz) stability estimates of the nonlinear seismic inverse problem.

Novel iterative methods, conditional Lipschitz stability and convergence

Forward modelling and optimization to match waveforms is an approach that was introduced in seismology more than four decades ago. Even though this approach, commonly using implementations based on the ‘adjoint state’ method, is being rather widely applied today, no results on the recovery of coefficients of models, however, were obtained. The first fundamental question to be addressed is whether the inverse boundary value problem corresponding with idealized seismic data has a unique solution. The second question is whether explicit reconstructions can be obtained, and the third question is what one can say about the case of partial data, also subject to sampling. The fourth question concerns the incorporation of inaccurate data and uncertainty quantification. All these questions are intimately interconnected. The second question typically involves iterative methods (‘full waveform inversion’) and leads to a local result.

We succeeded in obtaining conditional convergence of iterative methods results; we obtained explicit expressions for the radii of convergence [96]. The convergence relies on a Hölder, or Lipschitz, stability estimate for (essentially well-posedness of) the inverse problem. These insights have a significant impact on our understanding and develop- ment of optimization techniques very different from optimization with traditional regularization strategies.

  • Lipschitz stability estimates: scalar waves. The time-harmonic or ‘frequency-domain’ formulation lends itself to obtain conditional Lipschitz stability estimates for the nonlinear seismic inverse boundary value problem, even assuming partial boundary data [103]. The data are the Dirichlet-to-Neumann map. Uniqueness of recovery for this inverse problem in the case of scalar waves has been known since the late 1980s under minimal regularity assumptions. We established these estimates in the case of piecewise constant models with underlying domain partitionings containing discontinuities.

    The Lipschitz stability constant grows exponentially with the number of subdomains in the domain partitioning. The stability constant appears in the expressions for the radii of convergence, in a way such that the radius tends to zero as this constant becomes large. We obtained a sharp estimate of the stability constant in terms of the number of subdomains in the domain partitioning [preprint 7]. This estimate plays a critical role in the development of a multi-level (multi-scale) scheme.

    Moreover, we obtained the Lipschitz stable recovery of polyhedral interfaces and – if the domain partitioning in the above is an (unstructured) tetrahedral mesh – the Lipschitz stable recovery of the number of elements in and the shape of the mesh [preprint 8]. Thus we fundamentally relaxed knowledge about the parametrization.

    The Dirichlet-to-Neumann map is related to a single layer potential operator, the kernel of which is directly related to seismic data.

  • Time-harmonic elastic waves. We obtained the unique recovery of Lame´ parameters and density that are piece- wise constant on an unstructured tetrahedral mesh from the Neumann-to-Dirichlet map (which corresponds with vibroseis data). We also obtained Lipschitz stability.
  • Multi-level, multi-scale schemes. Based on the conditional stability estimates we introduce convex subsets, fol- lowing a refinement of domain partitionings. We designed an iterative method which converges to the point in each subset closest to the unique solution. The distance, or approximation error, depends on frequency, and an estimate can be obtained. Mitigating the growth of the stability constants with the approximation errors controlled by frequency, we obtained a multi-level (nonlinear projected steepest descent iteration) scheme with precise conditions for convergence [109] defining a nonlinear hierarchical compressed reconstruction. Through the compression and corresponding approximation we incorporate inaccuracies in the data. Other more ad- vanced methods, in particular the Gauss-Newton method, can also be adapted for this purpose.
  • Attenuation. The Lipschitz stability estimates hold for complex (frequency dependent) coefficients and, hence, the time-harmonic formulation allows in principle to obtain localized information about the (reciprocal) Q factor, and spatial variability in physical dispersion. This might open the way to the further studying of attenuation in localized regions in the mantle, such as in the asthenosphere.

Hyperbolic inverse boundary value problem

The time-domain formulation of the inverse boundary value problem for scalar and elastic waves involves different techniques and strategies. Naturally, in this formulation, one can introduce time windows identifying structure in and emphasizing parts of the data. We have been studying this problem for coefficients with structure. On the one hand, in the case of piecewise smooth coefficients with interfaces of unknown shapes, we use connections with geometric inverse problems with corresponding conditions; on the other hand, we establish explicit connections with the time- harmonic inverse boundary value problem, a multi-level strategy, and conditional high-frequency stability estimates.

Conditional Lipschitz Stability, Sparsity, Convergence Radius

Nonlinear inverse (boundary value) problems: Conditional well-posedness, iterative methods

Part of my research program has been devoted to developing new, comprehensive quantitative probes of our planet’s interior which include novel convergent nonlinear iterative techniques (surpassing so-called standard ‘adjoint state’ methods). One of our most recent breakthroughs concerns the understanding of the fundamental role of conditional Hölder (or Lipschitz) stability estimates of the nonlinear seismic inverse problem.

Novel iterative methods, conditional Lipschitz stability and convergence

Forward modelling and optimization to match waveforms is an approach that was introduced in seismology more than four decades ago. Even though this approach, commonly using implementations based on the ‘adjoint state’ method, is being rather widely applied today, no results on the recovery of coefficients of models, however, were obtained. The first fundamental question to be addressed is whether the inverse boundary value problem corresponding with idealized seismic data has a unique solution. The second question is whether explicit reconstructions can be obtained, and the third question is what one can say about the case of partial data, also subject to sampling. The fourth question concerns the incorporation of inaccurate data and uncertainty quantification. All these questions are intimately interconnected. The second question typically involves iterative methods (‘full waveform inversion’) and leads to a local result.

We succeeded in obtaining conditional convergence of iterative methods results; we obtained explicit expressions for the radii of convergence [96]. The convergence relies on a Hölder, or Lipschitz, stability estimate for (essentially well-posedness of) the inverse problem. These insights have a significant impact on our understanding and develop- ment of optimization techniques very different from optimization with traditional regularization strategies.

  • Lipschitz stability estimates: scalar waves. The time-harmonic or ‘frequency-domain’ formulation lends itself to obtain conditional Lipschitz stability estimates for the nonlinear seismic inverse boundary value problem, even assuming partial boundary data [103]. The data are the Dirichlet-to-Neumann map. Uniqueness of recovery for this inverse problem in the case of scalar waves has been known since the late 1980s under minimal regularity assumptions. We established these estimates in the case of piecewise constant models with underlying domain partitionings containing discontinuities.The Lipschitz stability constant grows exponentially with the number of subdomains in the domain partitioning. The stability constant appears in the expressions for the radii of convergence, in a way such that the radius tends to zero as this constant becomes large. We obtained a sharp estimate of the stability constant in terms of the number of subdomains in the domain partitioning [preprint 7]. This estimate plays a critical role in the development of a multi-level (multi-scale) scheme.

    Moreover, we obtained the Lipschitz stable recovery of polyhedral interfaces and – if the domain partitioning in the above is an (unstructured) tetrahedral mesh – the Lipschitz stable recovery of the number of elements in and the shape of the mesh [preprint 8]. Thus we fundamentally relaxed knowledge about the parametrization.

    The Dirichlet-to-Neumann map is related to a single layer potential operator, the kernel of which is directly related to seismic data.

  • Time-harmonic elastic waves. We obtained the unique recovery of Lame´ parameters and density that are piece- wise constant on an unstructured tetrahedral mesh from the Neumann-to-Dirichlet map (which corresponds with vibroseis data). We also obtained Lipschitz stability.
  • Multi-level, multi-scale schemes. Based on the conditional stability estimates we introduce convex subsets, fol- lowing a refinement of domain partitionings. We designed an iterative method which converges to the point in each subset closest to the unique solution. The distance, or approximation error, depends on frequency, and an estimate can be obtained. Mitigating the growth of the stability constants with the approximation errors controlled by frequency, we obtained a multi-level (nonlinear projected steepest descent iteration) scheme with precise conditions for convergence [109] defining a nonlinear hierarchical compressed reconstruction. Through the compression and corresponding approximation we incorporate inaccuracies in the data. Other more ad- vanced methods, in particular the Gauss-Newton method, can also be adapted for this purpose.
  • Attenuation. The Lipschitz stability estimates hold for complex (frequency dependent) coefficients and, hence, the time-harmonic formulation allows in principle to obtain localized information about the (reciprocal) Q factor, and spatial variability in physical dispersion. This might open the way to the further studying of attenuation in localized regions in the mantle, such as in the asthenosphere.

Hyperbolic inverse boundary value problem

The time-domain formulation of the inverse boundary value problem for scalar and elastic waves involves different techniques and strategies. Naturally, in this formulation, one can introduce time windows identifying structure in and emphasizing parts of the data. We have been studying this problem for coefficients with structure. On the one hand, in the case of piecewise smooth coefficients with interfaces of unknown shapes, we use connections with geometric inverse problems with corresponding conditions; on the other hand, we establish explicit connections with the time- harmonic inverse boundary value problem, a multi-level strategy, and conditional high-frequency stability estimates.

Earth as an unstructured mesh

Nonlinear inverse (boundary value) problems: Conditional well-posedness, iterative methods

Part of my research program has been devoted to developing new, comprehensive quantitative probes of our planet’s interior which include novel convergent nonlinear iterative techniques (surpassing so-called standard ‘adjoint state’ methods). One of our most recent breakthroughs concerns the understanding of the fundamental role of conditional Hölder (or Lipschitz) stability estimates of the nonlinear seismic inverse problem.

Novel iterative methods, conditional Lipschitz stability and convergence

Forward modelling and optimization to match waveforms is an approach that was introduced in seismology more than four decades ago. Even though this approach, commonly using implementations based on the ‘adjoint state’ method, is being rather widely applied today, no results on the recovery of coefficients of models, however, were obtained. The first fundamental question to be addressed is whether the inverse boundary value problem corresponding with idealized seismic data has a unique solution. The second question is whether explicit reconstructions can be obtained, and the third question is what one can say about the case of partial data, also subject to sampling. The fourth question concerns the incorporation of inaccurate data and uncertainty quantification. All these questions are intimately interconnected. The second question typically involves iterative methods (‘full waveform inversion’) and leads to a local result.

We succeeded in obtaining conditional convergence of iterative methods results; we obtained explicit expressions for the radii of convergence [96]. The convergence relies on a Hölder, or Lipschitz, stability estimate for (essentially well-posedness of) the inverse problem. These insights have a significant impact on our understanding and develop- ment of optimization techniques very different from optimization with traditional regularization strategies.

  • Lipschitz stability estimates: scalar waves. The time-harmonic or ‘frequency-domain’ formulation lends itself to obtain conditional Lipschitz stability estimates for the nonlinear seismic inverse boundary value problem, even assuming partial boundary data [103]. The data are the Dirichlet-to-Neumann map. Uniqueness of recovery for this inverse problem in the case of scalar waves has been known since the late 1980s under minimal regularity assumptions. We established these estimates in the case of piecewise constant models with underlying domain partitionings containing discontinuities.

    The Lipschitz stability constant grows exponentially with the number of subdomains in the domain partitioning. The stability constant appears in the expressions for the radii of convergence, in a way such that the radius tends to zero as this constant becomes large. We obtained a sharp estimate of the stability constant in terms of the number of subdomains in the domain partitioning [preprint 7]. This estimate plays a critical role in the development of a multi-level (multi-scale) scheme.

    Moreover, we obtained the Lipschitz stable recovery of polyhedral interfaces and – if the domain partitioning in the above is an (unstructured) tetrahedral mesh – the Lipschitz stable recovery of the number of elements in and the shape of the mesh [preprint 8]. Thus we fundamentally relaxed knowledge about the parametrization.

    The Dirichlet-to-Neumann map is related to a single layer potential operator, the kernel of which is directly related to seismic data.

  • Time-harmonic elastic waves. We obtained the unique recovery of Lame´ parameters and density that are piece- wise constant on an unstructured tetrahedral mesh from the Neumann-to-Dirichlet map (which corresponds with vibroseis data). We also obtained Lipschitz stability.
  • Multi-level, multi-scale schemes. Based on the conditional stability estimates we introduce convex subsets, fol- lowing a refinement of domain partitionings. We designed an iterative method which converges to the point in each subset closest to the unique solution. The distance, or approximation error, depends on frequency, and an estimate can be obtained. Mitigating the growth of the stability constants with the approximation errors controlled by frequency, we obtained a multi-level (nonlinear projected steepest descent iteration) scheme with precise conditions for convergence [109] defining a nonlinear hierarchical compressed reconstruction. Through the compression and corresponding approximation we incorporate inaccuracies in the data. Other more ad- vanced methods, in particular the Gauss-Newton method, can also be adapted for this purpose.
  • Attenuation. The Lipschitz stability estimates hold for complex (frequency dependent) coefficients and, hence, the time-harmonic formulation allows in principle to obtain localized information about the (reciprocal) Q factor, and spatial variability in physical dispersion. This might open the way to the further studying of attenuation in localized regions in the mantle, such as in the asthenosphere.

Hyperbolic inverse boundary value problem

The time-domain formulation of the inverse boundary value problem for scalar and elastic waves involves different techniques and strategies. Naturally, in this formulation, one can introduce time windows identifying structure in and emphasizing parts of the data. We have been studying this problem for coefficients with structure. On the one hand, in the case of piecewise smooth coefficients with interfaces of unknown shapes, we use connections with geometric inverse problems with corresponding conditions; on the other hand, we establish explicit connections with the time- harmonic inverse boundary value problem, a multi-level strategy, and conditional high-frequency stability estimates.

Wave Propagation and Scattering in Anisotropic Solids

Fast algorithms, massively parallel computing

Direct structured solvers and time-harmonic waves

In the inverse 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.

Segmentation Driven Domain Partitioning

Nonlinear inverse (boundary value) problems: Conditional well-posedness, iterative methods

Part of my research program has been devoted to developing new, comprehensive quantitative probes of our planet’s interior which include novel convergent nonlinear iterative techniques (surpassing so-called standard ‘adjoint state’ methods). One of our most recent breakthroughs concerns the understanding of the fundamental role of conditional Hölder (or Lipschitz) stability estimates of the nonlinear seismic inverse problem.

Novel iterative methods, conditional Lipschitz stability and convergence

Forward modelling and optimization to match waveforms is an approach that was introduced in seismology more than four decades ago. Even though this approach, commonly using implementations based on the ‘adjoint state’ method, is being rather widely applied today, no results on the recovery of coefficients of models, however, were obtained. The first fundamental question to be addressed is whether the inverse boundary value problem corresponding with idealized seismic data has a unique solution. The second question is whether explicit reconstructions can be obtained, and the third question is what one can say about the case of partial data, also subject to sampling. The fourth question concerns the incorporation of inaccurate data and uncertainty quantification. All these questions are intimately interconnected. The second question typically involves iterative methods (‘full waveform inversion’) and leads to a local result.

We succeeded in obtaining conditional convergence of iterative methods results; we obtained explicit expressions for the radii of convergence [96]. The convergence relies on a Hölder, or Lipschitz, stability estimate for (essentially well-posedness of) the inverse problem. These insights have a significant impact on our understanding and develop- ment of optimization techniques very different from optimization with traditional regularization strategies.

  • Lipschitz stability estimates: scalar waves. The time-harmonic or ‘frequency-domain’ formulation lends itself to obtain conditional Lipschitz stability estimates for the nonlinear seismic inverse boundary value problem, even assuming partial boundary data [103]. The data are the Dirichlet-to-Neumann map. Uniqueness of recovery for this inverse problem in the case of scalar waves has been known since the late 1980s under minimal regularity assumptions. We established these estimates in the case of piecewise constant models with underlying domain partitionings containing discontinuities.

    The Lipschitz stability constant grows exponentially with the number of subdomains in the domain partitioning. The stability constant appears in the expressions for the radii of convergence, in a way such that the radius tends to zero as this constant becomes large. We obtained a sharp estimate of the stability constant in terms of the number of subdomains in the domain partitioning [preprint 7]. This estimate plays a critical role in the development of a multi-level (multi-scale) scheme.

    Moreover, we obtained the Lipschitz stable recovery of polyhedral interfaces and – if the domain partitioning in the above is an (unstructured) tetrahedral mesh – the Lipschitz stable recovery of the number of elements in and the shape of the mesh [preprint 8]. Thus we fundamentally relaxed knowledge about the parametrization.

    The Dirichlet-to-Neumann map is related to a single layer potential operator, the kernel of which is directly related to seismic data.

  • Time-harmonic elastic waves. We obtained the unique recovery of Lame´ parameters and density that are piece- wise constant on an unstructured tetrahedral mesh from the Neumann-to-Dirichlet map (which corresponds with vibroseis data). We also obtained Lipschitz stability.
  • Multi-level, multi-scale schemes. Based on the conditional stability estimates we introduce convex subsets, fol- lowing a refinement of domain partitionings. We designed an iterative method which converges to the point in each subset closest to the unique solution. The distance, or approximation error, depends on frequency, and an estimate can be obtained. Mitigating the growth of the stability constants with the approximation errors controlled by frequency, we obtained a multi-level (nonlinear projected steepest descent iteration) scheme with precise conditions for convergence [109] defining a nonlinear hierarchical compressed reconstruction. Through the compression and corresponding approximation we incorporate inaccuracies in the data. Other more ad- vanced methods, in particular the Gauss-Newton method, can also be adapted for this purpose.
  • Attenuation. The Lipschitz stability estimates hold for complex (frequency dependent) coefficients and, hence, the time-harmonic formulation allows in principle to obtain localized information about the (reciprocal) Q factor, and spatial variability in physical dispersion. This might open the way to the further studying of attenuation in localized regions in the mantle, such as in the asthenosphere.

Hyperbolic inverse boundary value problem

The time-domain formulation of the inverse boundary value problem for scalar and elastic waves involves different techniques and strategies. Naturally, in this formulation, one can introduce time windows identifying structure in and emphasizing parts of the data. We have been studying this problem for coefficients with structure. On the one hand, in the case of piecewise smooth coefficients with interfaces of unknown shapes, we use connections with geometric inverse problems with corresponding conditions; on the other hand, we establish explicit connections with the time- harmonic inverse boundary value problem, a multi-level strategy, and conditional high-frequency stability estimates.