A preconditioned 3-D multi-region fast multipole solver for seismic wave propagation in complex geometries

International audienceThe analysis of seismic wave propagation and amplification in complex geological structures requires efficient numerical methods. In this article, following up on recent studies devoted to the formulation, implementation and evaluation of 3-D single- and multi-region elastodynamic fast multipole boundary element methods (FM-BEMs), a simple preconditioning strategy is proposed. Its efficiency is demonstrated on both the single- and multi-region versions using benchmark examples (scattering of plane waves by canyons and basins). Finally, the preconditioned FM-BEM is applied to the scattering of plane seismic waves in an actual configuration (alpine basin of Grenoble, France), for which the high velocity contrast is seen to significantly affect the overall efficiency of the multi-region FM-BEM

The method of polarized traces for the 2D Helmholtz equation

We present a solver for the 2D high-frequency Helmholtz equation in heterogeneous acoustic media, with online parallel complexity that scales optimally as O(NL), where N is the number of volume unknowns, and L is the number of processors, as long as L grows at most like a small fractional power of N. The solver decomposes the domain into layers, and uses transmission conditions in boundary integral form to explicitly define "polarized traces", i.e., up- and down-going waves sampled at interfaces. Local direct solvers are used in each layer to precompute traces of local Green's functions in an embarrassingly parallel way (the offline part), and incomplete Green's formulas are used to propagate interface data in a sweeping fashion, as a preconditioner inside a GMRES loop (the online part). Adaptive low-rank partitioning of the integral kernels is used to speed up their application to interface data. The method uses second-order finite differences. The complexity scalings are empirical but motivated by an analysis of ranks of off-diagonal blocks of oscillatory integrals. They continue to hold in the context of standard geophysical community models such as BP and Marmousi 2, where convergence occurs in 5 to 10 GMRES iterations. While the parallelism in this paper stems from decomposing the domain, we do not explore the alternative of parallelizing the systems solves with distributed linear algebra routines. Keywords: Domain decomposition; Helmholtz equation; Integral equations; High-frequency; Fast methodsUnited States. Air Force Office of Scientific Research (Grant FA9550-15-1-0078)United States. Office of Naval Research (Grant N00014-13-1-0403)National Science Foundation (U.S.) (Grant DMS-1255203

Efficient discretisation and domain decomposition preconditioners for incompressible fluid mechanics

Solving the linear elasticity and Stokes equations by an optimal domain decomposition method derived algebraically involves the use of non standard interface conditions whose discretisation is not trivial. For this reason the use of approximation methods such as hybrid discontinuous Galerkin appears as an appropriate strategy: on the one hand they provide the best compromise in terms of the number of degrees of freedom in between standard continuous and discontinuous Galerkin methods, and on the other hand the degrees of freedom used in the non standard interface conditions are naturally defined at the boundary between elements. In this manuscript we present the coupling between a well chosen discretisation method (hybrid discontinuous Galerkin) and a novel and efficient domain decomposition method to solve the Stokes system. An analysis of the boundary value problem with non standard condition is provided as well as the numerical evidence showing the advantages of the new method. Furthermore, we present and analyse a stabilisation method for the presented discretisation that allows the use of the same polynomial degrees for velocity and pressure discrete spaces. The original definition of the domain decomposition preconditioners is one-level, this is, the preconditioner is built only using the solution of local problems. This has the undesired consequence that the results are not scalable, it means that the number of iterations needed to reach convergence increases with the number of subdomains. This is the reason why we have also introduced, and tested numerically, two-level preconditioners. Such preconditioners use a coarse space in their construction. We consider two finite element discretisations, namely, the hybrid discontinuous Galerkin and Taylor-Hood discretisations for the nearly incompressible elasticity problems and Stokes equations.Solving the linear elasticity and Stokes equations by an optimal domain decomposition method derived algebraically involves the use of non standard interface conditions whose discretisation is not trivial. For this reason the use of approximation methods such as hybrid discontinuous Galerkin appears as an appropriate strategy: on the one hand they provide the best compromise in terms of the number of degrees of freedom in between standard continuous and discontinuous Galerkin methods, and on the other hand the degrees of freedom used in the non standard interface conditions are naturally defined at the boundary between elements. In this manuscript we present the coupling between a well chosen discretisation method (hybrid discontinuous Galerkin) and a novel and efficient domain decomposition method to solve the Stokes system. An analysis of the boundary value problem with non standard condition is provided as well as the numerical evidence showing the advantages of the new method. Furthermore, we present and analyse a stabilisation method for the presented discretisation that allows the use of the same polynomial degrees for velocity and pressure discrete spaces. The original definition of the domain decomposition preconditioners is one-level, this is, the preconditioner is built only using the solution of local problems. This has the undesired consequence that the results are not scalable, it means that the number of iterations needed to reach convergence increases with the number of subdomains. This is the reason why we have also introduced, and tested numerically, two-level preconditioners. Such preconditioners use a coarse space in their construction. We consider two finite element discretisations, namely, the hybrid discontinuous Galerkin and Taylor-Hood discretisations for the nearly incompressible elasticity problems and Stokes equations

Prediction of Acoustic Resonances in Core Volumes

In numerous industrial systems, acoustic resonances are synonymous with loud noise and dramatic structural vibrations. They take their origin from trapped weak pressure waves, but can sometimes develop into complex instabilities. With specific emphasis on turbomachines applications, two computational methods have been developed for the prediction of acoustic resonances in core volumes. The first method, called the averaged response function, consists of combined time-domain and frequency-domain approaches. A Favre-averaged RANS solver is used to model the response of a system to a controlled chirp excitation. The response is then analysed using well-known modal analysis tools in order to extract the system’s acoustic characteristics. The second method, called the Arnoldi method, focuses on the stability of the CFD solver. The system is transformed, using a linear Euler solver, into a classic matrix stability problem. The eigenpairs of the matrix, corresponding to the acoustic modeshapes of the system, are then extracted thanks to the iterative Arnoldi method. The two methods are validated and compared on a wide range of cases such as enclosures, open geometries and flow applications. Such a thorough study first provides the reader with a deeper insight into acoustic phenomena by considering classic acoustic examples such as the end correction concept, the Doppler effect and the ”lock-in“ phenomenon. This study also investigates the limitations and qualities of the implemented methods which are seen to behave very well when compared to theory and experiments. They give accurate results in predicting the three components of the acoustic resonance: the frequency, the damping and the modeshape. As a result, the methods implemented are considered to be mature and can be used to study either the complete acoustic map of the system across a wide range of frequencies or specific acoustic instabilities in a narrow frequency range

Application of the wavelet transform for sparse matrix systems and PDEs.

Thesis (M.Sc.)-University of KwaZulu-Natal, Westville, 2009.We consider the application of the wavelet transform for solving sparse matrix systems and partial differential equations. The first part is devoted to the theory and algorithms of wavelets. The second part is concerned with the sparse representation of matrices and well-known operators. The third part is directed to the application of wavelets to partial differential equations, and to sparse linear systems resulting from differential equations. We present several numerical examples and simulations for the above cases

Microwave Imaging of Brain Stroke:Contributions to Modeling and Inverse Problem Resolution

Brain stroke is an age-related illness which has become a major issue in our ageing societies. Early diagnosis and treatment are of high importance for the full recovery of the patient, as reminded in Anglo-Saxon countries by the abbreviation FAST (Face, Arm, Speech, Time) referring to both the four major visible signs and the necessity to act fast. In this respect, Computed Tomography (CT) and Nuclear Magnetic Resonance (NMR) imaging are key diagnostic tools in clinical practice. Unfortunately, not only these modalities can neither be transported nor rapidly usable, which would allow early treatment (especially in rural environments), but also cannot be brought to the bedside of the patient to monitor the evolution of the disease. Microwave Imaging (MWI) is a potential candidate to provide fast and accurate diagnostic insights for brain stroke pathological states. The head of the patient is illuminated with low-power microwave waveforms (non-ionizing radiations), whose backscattered signals are used to generate either images of its internal structures, distributions, patterns and shapes (qualitative imaging) or directly its physical parameters such as the dielectric contrast and the permittivity values (quantitative imaging). The technology relies on the high sensitivity of microwaves on the water content of tissues to allow for the discrimination between pathological and healthy regions. This thesis focuses on both the forward modeling of the electromagnetic phenomena arising in biological tissues and the inverse scattering problem for imaging in the differential MWI (dMWI) scenario for brain stroke monitoring. It is intrinsically interdisciplinary as it requires knowledge in Biology, Medicine, Physics, Chemistry, and Engineering. In order to investigate the challenges arising in brain MWI, it is crucial to have accurate and efficient solvers to model electromagnetic (EM) fields at UHF/SHF-bands. The head is a distributed, heterogeneous, and lossy scatterer for which existing solvers are known to struggle at higher frequencies. Volume Integral Equation (VIE) formulations and MultiGrid (MG) approaches are investigated to find the actual solution of the field distributions for large scale problems. The EM modeling also permits to analyze the feasibility of brain MWI, which depends on the power transmission from the antennas towards the human brain. In order to estimate this transmission, simplified but still representative models, including intermediate layers -skin, fat, bone, and CerebroSpinal Fluid (CSF) - of the head, are proposed in the framework of simulations (analytical tools) and experimental validations (3D printed head phantom). For the imaging task, the physics of the EM scattering, leads to complex non-linear inverse scattering problems (consisting in retrieving from a set of field measurements the physical parameters which produced them) for which reliable assumptions and approximations must be found. For brain MWI, estimating and quantifying the degree of non-linearity allows for determining the scope of application of existing algorithms, for which different regularizers are applied. Modeling and inverse problem resolution for brain MWI investigated in the present work are ultimately meant to contribute to the development of a technology dedicated to brain stroke detection, differentiation, and monitoring

Méthodes quasi-optimales pour la résolution des équations intégrales de frontière en électromagnétisme

Il existe une grande quantité de méthodes numériques adaptées d’une part à la modélisation, et d'autre part à la résolution des équations de Maxwell. En particulier, la méthode des éléments nis de frontière (BEM), ou méthode des Moments (MoM), semble appropriée pour la mise en équation des phénomènes de diffraction par des objets parfaitement conducteurs, en limitant le cadre de l'étude à la frontière entre l'objet diffractant et le milieu extérieur. Cette méthode mène systématiquement à la résolution d’un système linéaire dense, que nous parvenons à compresser en l'approchant numériquement par une matrice hiérarchique creuse, appelée H-matrice. Cette approximation peut être complétée d'une ré-agglomération permettant d'améliorer la sparsité de la H-matrice et ainsi d'optimiser davantage la résolution du système traité. La hiérarchisation du système s'effectue en considérant la matrice traitée par blocs, que l'on peut ou non compresser selon une condition d'admissibilité. L'Approximation en Croix Adaptative (ACA) ou l'Approximation en Croix Hybride (HCA) sont deux méthodes de compression que l'on peut alors appliquer aux blocs admissibles. Il existe une grande quantité de méthodes numériques adaptées d’une part à la modélisation, et d'autre part à la résolution des équations de Maxwell. En particulier, la méthode des éléments finis de frontière (BEM), ou méthode des Moments (MoM), semble appropriée pour la mise en équation des phénomènes de diffraction par des objets parfaitement conducteurs, en limitant le cadre de l'étude à la frontière entre l'objet diffractant et le milieu extérieur. Cette méthode mène systématiquement à la résolution d’un système linéaire dense, que nous parvenons à compresser en l'approchant numériquement par une matrice hiérarchique creuse, appelée H-matrice. Cette approximation peut être complétée d'une ré-agglomération permettant d'améliorer la sparsité de la H-matrice et ainsi d'optimiser davantage la résolution du système traité. La hiérarchisation du système s'effectue en considérant la matrice traitée par blocs, que l'on peut ou non compresser selon une condition d'admissibilité. L'Approximation en Croix Adaptative (ACA) ou l'Approximation en Croix Hybride (HCA) sont deux méthodes de compression que l'on peut alors appliquer aux blocs admissibles. Le travail de cette thèse consiste dans un premier temps à valider le format H-matrice en 2D et en 3D en utilisant l'ACA, puis d'y appliquer la méthode HCA, encore peu exploitée. Nous pouvons alors résoudre le système linéaire issu de la BEM en utilisant différents solveurs, directs ou non, adaptés au format hiérarchique. En particulier, nous pourrons constater l'efficacité du préconditionnement LU hiérarchique sur un solveur itératif. Nous pourrons alors appliquer ce formalisme au cas des surfaces rugueuses ou encore des fibres à cristaux photoniques (PCF). Il sera également possible de paralléliser certaines opérations sur architecture partagée afin de réduire de nouveau le coût temporel de la résolution