Search

International audience; Given a full-rank matrix $A \in \mathbb{R}^{m\times n}$ ($m\geq n$), we consider a special class of linear systems $A^T\! Ax=A^T\! b+c$ with $x, c \in \mathbb{R}^{n}$ and $b \in \mathbb{R}^{m}$, which we refer to as the extended normal equations. The occurrence of $c$ gives rise to a problem with a different conditioning from the standard normal equations and prevents direct application of standard methods for least squares. Hence, we seek more insights on theoretical and practical aspects of the solution of such problems. We propose an explicit formula for the structured condition number, which allows us to compute a more accurate estimate of the forward error than the standard one used for generic linear systems, which does not take into account the structure of the perturbations. The relevance of our estimate is shown on a set of synthetic test problems. Then, we propose a new iterative solution method that, as in the case of normal equations, takes advantage of the structure of the system to avoid unstable computations such as forming $A^T\! A$ explicitly. Numerical experiments highlight the increased robustness and accuracy of the proposed method compared to standard iterative methods. It is also found that the new method can compare to standard direct methods in terms of solution accuracy.

SIAM: Society for Industrial and Applied Mathematics; International audience; We propose a new family of multilevel methods for unconstrained minimization. The resulting strategies are multilevel extensions of high-order optimization methods based on q-order Taylor models (with q ≥ 1) that have been recently proposed in the literature. The use of high-order models, while decreasing the worst-case complexity bound, makes these methods computationally more expensive. Hence, to counteract this effect, we propose a multilevel strategy that exploits a hierarchy of problems of decreasing dimension, still approximating the original one, to reduce the global cost of the step computation. A theoretical analysis of the family of methods is proposed. Specifically, local and global convergence results are proved and a complexity bound to reach first order stationary points is also derived. A multilevel version of the well known adaptive method based on cubic regularization (ARC, corresponding to q = 2 in our setting) has been implemented. Numerical experiments clearly highlight the relevance of the new multilevel approach leading to considerable computational savings in terms of floating point operations compared to the classical one-level strategy.

We consider the design of structure-preserving discretization methods for the solution of systems of boundary controlled Partial Differential Equations (PDEs) thanks to the port-Hamiltonian formalism. We first provide a novel general structure of infinite-dimensional port-Hamiltonian systems (pHs) for which the Partitioned Finite Element Method (PFEM) straightforwardly applies. The proposed strategy is applied to abstract multidimensional linear hyperbolic and parabolic systems of PDEs. Then we show that instructional model problems based on the wave equation, Mindlin equation and heat equation fit within this unified framework. Secondly we introduce the ongoing project SCRIMP (Simulation and ContRol of Interactions in Multi-Physics) developed for the numerical simulation of infinite-dimensional pHs. SCRIMP notably relies on the FEniCS open-source computing platform for the finite element spatial discretization. Finally, we illustrate how to solve the considered model problems within this framework by carefully explaining the methodology. As additional support, companion interactive Jupyter notebooks are available., 33 pages

International audience; Many applications in structural mechanics require the numerical solution of sequences of linear systems typically issued from a finite element discretization of the governing equations on fine meshes. The method of Lagrange multipliers is often used to take into account mechanical constraints. The resulting matrices then exhibit a saddle point structure and the iterative solution of such preconditioned linear systems is considered as challenging. A popular strategy is then to combine preconditioning and deflation to yield an efficient method.We propose an alternative that is applicable to the general case and not only to matrices with a saddle point structure. In this approach, we consider to update an existing algebraic or application-based preconditioner, using specific available information exploiting the knowledge of an approximate invariant subspace or of matrix-vector products. The resulting preconditioner has the form of a limited memory quasi-Newton matrix and requires a small number of linearly independent vectors. Numerical experiments performed on three large-scale applications in elasticity highlight the relevance of the new approach. We show that the proposed method outperforms the deflation method when considering sequences of linear systems with varying matrices.

International audience; In this paper, we propose a new multilevel Levenberg–Marquardt optimizer for the training of artificial neural networks with quadratic loss function. This setting allows us to get further insight into the potential of multilevel optimization methods. Indeed, when the least squares problem arises from the training of artificial neural networks,the variables subject to optimization are not related by any geometrical constraints and the standard interpolation and restriction operators cannot be employed any longer. A heuristic, inspired by algebraic multigrid methods, is then proposed to construct the multilevel transfer operators. We test the new optimizer on an important application: the approximate solution of partial differential equations by means of artificial neural networks. The learning problem is formulated as a least squares problem, choosing the nonlinear residual of the equation as a loss function, whereas the multilevel method is employed as a training method. Numerical experiments show encouraging results related to the efficiency of the new multilevel optimization method compared to the corresponding one-level procedure in this context.

International audience; We propose a method for representing spatially correlated observation errors in variational data assimilation. The method is based on the numerical solution of a diffusion equation, a technique commonly used for representing spatially correlated background errors. The discretization of the pseudo‐time derivative of the diffusion equation is done implicitly using a backward Euler scheme. The solution of the resulting elliptic equation can be interpreted as a correlation operator whose kernel is a correlation function from the Matérn family.In order to account for the possibly heterogeneous distribution of observations, a spatial discretization technique based on the finite element method (FEM) is chosen where the observation locations are used to define the nodes of an unstructured mesh on which the diffusion equation is solved. By construction, the method leads to a convenient operator for the inverse of the observation‐error correlation matrix, which is an important requirement when applying it with standard minimization algorithms in variational data assimilation. Previous studies have shown that spatially correlated observation errors can also be accounted for by assimilating the observations together with their directional derivatives up to arbitrary order. In the continuous framework, we show that the two approaches are formally equivalent for certain parameter specifications. The FEM provides an appropriate framework for evaluating the derivatives numerically, especially when the observations are heterogeneously distributed.Numerical experiments are performed using a realistic data distribution from the Spinning Enhanced Visible and InfraRed Imager (SEVIRI). Correlations obtained with the FEM‐discretized diffusion operator are compared with those obtained using the analytical Matérn correlation model. The method is shown to produce an accurate representation of the target Matérn function in regions where the data are densely distributed. The presence of large gaps in the data distribution degrades the quality of the mesh and leads to numerical errors in the representation of the Matérn function. Strategies to improve the accuracy of the method in the presence of such gaps are discussed.

International audience; This paper is concerned with the approximation of the solution of partial differential equations by means of artificial neural networks. Here a feedforward neural network is used to approximate the solution of the partial differential equation. The learning problem is formulated as a least squares problem, choosing the residual of the partial differential equation as a loss function, whereas a multilevel Levenberg-Marquardt method is employed as a training method. This setting allows us to get further insight into the potential of multilevel methods. Indeed, when the least squares problem arises from the training of artificial neural networks, the variables subject to optimization are not related by any geometrical constraints and the standard interpolation and restriction operators cannot be employed any longer. A heuristic, inspired by algebraic multigrid methods, is then proposed to construct the multilevel transfer operators. Numerical experiments show encouraging results related to the efficiency of the new multilevel optimization method for the training of artificial neural networks, compared to the standard corresponding one-level procedure.

The numerical solution of saddle point systems has received a lot of attention over the past few years in a wide variety of applications such as constrained optimization, computational fluid dynamics and optimal control, to name a few. In this paper, we focus on the saddle point formulation of a large-scale variational data assimilation problem, where the computations involving the constraint blocks are supposed to be much more expensive than those related to the (1, 1) block of the saddle point matrix. New low-rank limited memory preconditioners exploiting the particular structure of the problem are proposed and analysed theoretically. Numerical experiments performed within the Object-Oriented Prediction System are presented to highlight the relevance of the proposed preconditioners.

Direct Numerical Simulation of turbulent flows is a computationally demanding problem that requires efficient parallel algorithms. We investigate the applicability of the time-parallel Parareal algorithm to an instructional case study related to the simulation of the decay of homogeneous isotropic turbulence in three dimensions. We combine a Parareal variant based on explicit time integrators and spatial coarsening with the space-parallel Hybrid Navier–Stokes solver. We analyse the performance of this space–time parallel solver with respect to speedup and quality of the solution. The results are compared with reference data obtained with a classical explicit integration, using an error analysis which relies on the energetic content of the solution. We show that a single Parareal iteration is able to reproduce with high fidelity the main statistical quantities characterizing the turbulent flow field.

We develop and analyse Neumann–Neumann methods for hp finite‐element approximations of scalar elliptic problems on geometrically refined boundary layer meshes in three dimensions. These are meshes that are highly anisotropic where the aspect ratio typically grows exponentially with the polynomial degree. The condition number of our preconditioners is shown to be independent of the aspect ratio of the mesh and of potentially large jumps of the coefficients. In addition, it only grows polylogarithmically with the polynomial degree, as in the case of p approximations on shape‐regular meshes. This work generalizes our previous one on two‐dimensional problems in Toselli & Vasseur (2003a, submitted to Numerische Mathematik, 2003c to appear in Comput. Methods Appl. Mech. Engng.) and the estimates derived here can be employed to prove condition number bounds for certain types of FETI methods., IMA Journal of Numerical Analysis, 24 (1), ISSN:0272-4979, ISSN:1464-3642

International audience; The spectral analysis is a basic tool to characterise the behaviour of any convection scheme. By nature, the solution projected onto the Fourier basis enables to estimate the dissipation and the dispersion associated with the spatial discretisation of the hyperbolic linear problem. In this paper, we wish to revisit such analysis, focusing attention on two key points. The first point concerns the effects of time integration on the spectral analysis. It is shown with standard high-order Finite Difference schemes dedicated to aeroacoustics that the time integration has an effect on the required number of points per wavelength. The situation depends on the choice of the coupled schemes (one for time integration, one for space derivative and one for the filter) and here, the compact scheme with its eighth-order filter seems to have a better spectral accuracy than the considered dispersion relation preserving scheme with its associated filter, especially in terms of dissipation. Secondly, such a coupled space time approach is applied to the new class of high-order spectral discontinuous approaches, focusing especially on the Spectral Difference method. A new way to address the specific spectral behaviour of the scheme is introduced first for wavenumbers in [0,pi], following the Matrix Power method. For wavenumbers above pi, an aliasing phenomenon always occurs but it is possible to understand and to control the aliasing of the signal. It is shown that aliasing depends on the polynomial degree and on the number of time steps. A new way to define dissipation and dispersion is introduced and applied to wavenumbers larger than it. Since the new criteria recover the previous results for wavenumbers below it, the new proposed approach is an extension of all the previous ones dealing with dissipation and dispersion errors. At last, since the standard Finite Difference schemes can serve as reference solution for their capability in aeroacoustics, it is shown that the Spectral Difference method is as accurate as (or even more accurate) than the considered Finite Difference schemes.

International audience; We propose a ''constrained'' least-squares approach to estimate regional maps of equivalent-water heights by inverting GRACE-based potential anomalies at satellite altitude. According to the energy integral method, the anomalies of difference of geopotential between the two GRACE vehicles are derived from along-track K-Band Range-Rate (KBRR) residuals that correspond mainly to the continental water storage changes, once a priori known accelerations (i.e. static field, polar movements, atmosphere and ocean masses including tides) are removed during the orbit adjustment process. Newton's first law merely enables the Difference of Potential Anomalies from accurate KBRR data and the equivalent-water heights to be recovered. Spatial constraints versus spherical distance between elementary surface tiles are introduced to stabilize the linear system to cancel the effects of the north-south striping. Unlike the ''mascons'' approach, no basis of orthogonal functions (e.g., spherical harmonics) is used, so that the proposed regional method does not suffer from drawbacks related to any spectrum truncation. Time series of 10-day regional maps over South America for 2006-2009 also prove to be consistent with independent datasets, namely the outputs of hydrological models, ''mascons'' and global GRACE solutions.

We study methods for the numerical solution of the Helmholtz equation for two-dimensional applications in geophysics. The common framework of the iterative methods in our study is a combination of an inner iteration with a geometric multigrid method used as a preconditioner and an outer iteration with a Krylov subspace method. The preconditioning system is based on either a pure or shifted Helmholtz operator. A multigrid iteration is used to approximate the inverse of this operator. The proposed solution methods are evaluated on a complex benchmark in geophysics involving highly variable coefficients and high wavenumbers. We compare this preconditioned iterative method with a direct method and a hybrid method that combines our iterative approach with a direct method on a reduced problem. We see that the hybrid method outperforms both the iterative and the direct approach.

International audience; In this paper we propose a new way to compute a rough approximation solution, to be later used as a warm starting point in a more refined optimization process, for a challenging global optimization problem related to Earth imaging in geophysics. The warm start con- sists of a velocity model that approximately solves a full-waveform inverse problem at low frequency. Our motivation arises from the availability of massively parallel computing plat- forms and the natural parallelization of evolution strategies as global optimization methods for continuous variables.Our first contribution consists of developing a new and efficient parametrization of the velocity models to significantly reduce the dimension of the original optimization space. Our second contribution is to adapt a class of evolution strategies to the specificity of the physical problem at hands where the objective function evaluation is known to be the most expen- sive computational part. A third contribution is the development of a parallel evolution strategy solver, taking advantage of a recently proposed modification of these class of evolu- tionary methods that ensures convergence and promotes better performance under moderate budgets.The numerical results presented demonstrate the effectiveness of the algorithm on a realistic 3D full-waveform inverse problem in geophysics. The developed numerical approach allows us to successfully solve an acoustic full-waveform inversion problem at low frequencies on a reasonable number of cores of a distributed memory computer.

International audience; We propose a recursive Kalman filtering approach to map regional spatio-temporal variations of terrestrial water mass over large continental areas, such as South America. Instead of correcting hydrology model outputs by the GRACE observations using a Kalman filter estimation strategy, regional 2-by-2 degree water mass solutions are constructed by integration of daily potential differences deduced from GRACE K-band range rate (KBRR) measurements. Recovery of regional water mass anomaly averages obtained by accumulation of information of daily noise-free simulated GRACE data shows that convergence is relatively fast and yields accurate solutions. In the case of cumulating real GRACE KBRR data contaminated by observational noise, the sequential method of step-by-step integration provides estimates of water mass variation for the period 2004–2011 by considering a set of suitable a priori error uncertainty parameters to stabilize the inversion. Spatial and temporal averages of the Kalman filter solutions over river basin surfaces are consistent with the ones computed using global monthly/10-day GRACE solutions from official providers CSR, GFZ and JPL. They are also highly correlated to in situ records of river discharges (70–95 %), especially for the Obidos station where the total outflow of the Amazon River is measured. The sparse daily coverage of the GRACE satellite tracks limits the time resolution of the regional Kalman filter solutions, and thus the detection of short-term hydrological events.

In this paper, we present extensive numerical tests showing the performance and robustness of a Balancing Neumann-Neumann method for the solution of algebraic linear systems arising from hp finite element approximations of scalar elliptic problems on geometrically refined boundary layer meshes in three dimensions. The numerical results are in good agreement with the theoretical bound for the condition number of the preconditioned operator derived in [Toselli and Vasseur, IMA J. Numer. Anal. 24 (2004) 123–156]. They confirm that the condition numbers are independent of the aspect ratio of the mesh and of potentially large jumps of the coefficients. Good results are also obtained for certain singularly perturbed problems. The condition numbers only grow polylogarithmically with the polynomial degree, as in the case of p approximations on shape-regular meshes [Pavarino, RAIRO: Model. Math. Anal. Numer. 31 (1997) 471–493]. This paper follows [Toselli and Vasseur, Comput. Methods Appl. Mech. Engrg. 192 (2003) 4551–4579] on two dimensional problems.

A family of dual-primal finite element tearing and interconnecting (FETI) methods for edge element approximations in two dimensions is proposed and analyzed. The primal constraints here are averages over subdomain edges. It is shown that the condition number of the corresponding method is independent of the number of substructures and grows only polylogarithmically with the number of unknowns associated with individual substructures. The estimate is also independent of the jumps of both of the coefficients of the original problem. Numerical results validating our theoretical bounds are given.

We present extensive numerical tests showing the performance and robustness of certain Balancing Neumann–Neumann and one-level FETI methods for the solution of algebraic linear systems arising from hp finite element approximations of scalar elliptic problems on geometrically refined boundary layer meshes in two dimensions. The numerical results confirm the theoretical results derived in [A. Toselli, X. Vasseur, Technical Report 02-15, Seminar fur Angewandte Mathematik, ETH, Zurich, September 2002]: the condition numbers are independent of the aspect ratio of the mesh and of potentially large jumps of the coefficients. In addition, they only grow polylogarithmically with the polynomial degree, as in the case of p approximations on shape-regular meshes.

We propose an unconstrained approach to recover regional time-variations of surface mass anomalies using Level-1 Gravity Recovery and Climate Experiment (GRACE) orbit observations, for reaching spatial resolutions of a few hundreds of kilometers. Potential differences between the twin GRACE vehicles are determined along short satellite tracks using the energy integral method (i.e., integration of orbit parameters vs. time) in a quasi-inertial terrestrial reference frame. Potential differences residuals corresponding mainly to changes in continental hydrology are then obtained after removing the gravitational effects of the known geophysical phenomena that are mainly the static part of the Earth’s gravity field and time-varying contributions to gravity (Sun, Moon, planets, atmosphere, ocean, tides, variations of Earth’s rotation axis) through ad hoc models. Regional surface mass anomalies are restored from potential difference anomalies of 10 to 30-day orbits onto 1◦ continental grids by regularization techniques based on singular value decomposition. Error budget analysis has been made by considering the important effects of spectrum truncation, the time length of observation (or spatial coverage of the data to invert) and for different levels of noise.

International audience; This work is concerned with the development and study of a minimum residual norm subspace method based on the generalized conjugate residual method with inner orthogonalization (GCRO) method that allows flexible preconditioning and deflated restarting for the solution of non-symmetric or non-Hermitian linear systems. First we recall the main features of flexible generalized minimum residual with deflated restarting (FGMRES-DR), a recently proposed algorithm of the same family but based on the GMRES method. Next we introduce the new inner-outer subspace method named FGCRO-DR. A theoretical comparison of both algorithms is then made in the case of flexible preconditioning. It is proved that FGCRO-DR and FGMRES-DR are algebraically equivalent if a collinearity condition is satisfied. While being nearly as expensive as FGMRES-DR in terms of computational operations per cycle, FGCRO-DR offers the additional advantage to be suitable for the solution of sequences of slowly changing linear systems (where both the matrix and right-hand side can change) through subspace recycling. Numerical experiments on the solution of multidimensional elliptic partial differential equations show the efficiency of FGCRO-DR when solving sequences of linear systems.

This work discusses the computation of the steady, turbulent, two-dimensional incompressible viscous flow on structured cell-centered collocated grids. A rather new computational approach – the so-called fully coupled procedure with defect correction technique – is presented as an alternative both to classical decoupled approaches (SIMPLE [Int. J. Heat Mass Transfer 15 (1972) 1787–1806], PISO [J. Comput. Phys. 62(1) (1986) 40–65] and variants) and to weakly coupled solution methods of Vanka type [Fifth symposium on Turbulent Shear flows, 1985. p. 20:27–32, J. Comput. Phys. 65 (1986) 138–158]. Its main features will be highlighted and detailed. This strategy is evaluated on two demanding test cases (the simulation of the separated flow past an AS240-B airfoil at high incidence (19°, Re=2×10 6 ) and the simulation of the wake flow behind a two-dimensional hill ( Re =60 000 )), for which documented experimental data are available. Both robustness and computational efficiency of the new approach are shown.

The Flexible GMRES (FGMRES [1]) and the GMRES with deflated restarting (GMRES-DR [2]) methods are two algorithms derived from GMRES [3], that are considered as powerful when solving large non hermitian systems of linear equations. GMRES-DR is a variant of GMRES with an improved restarting technique that maintains in the Krylov subspace harmonic Ritz vector from the previous restart. In situations where the convergence of restarted GMRES is slow and where the matrix has few eigenvalues close to the origin, this technique has proved very efficient. The new method that we propose is the Flexible GMRES with deflated restarting (FGMRES-DR [6]), which combines the two above mentioned algorithms in order to yield better performance. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)

R��sum��: Le travail effectu�� vise �� am��liorer la robustesse des strat��gies algorithmiques destin��es �� la r��solution num��rique des ��quations de Navier-Stokes en variables vitesse-pression dans le cadre de la mod��lisation d'��coulements de fluide visqueux incompressible. Sont ainsi adopt��es les formulations d��coupl��e et fortement coupl��e. L'objectif majeur de cette ��tude revient �� d��terminer, ��valuer et comparer diff��rentes techniques d'acc��l��ration de convergence au sein des formulations retenues. Acc��l��rer la convergence de toute m��thode d��coupl��e passe par l'emploi d'un solveur lin��aire robuste pour r��soudre l'��quation aux d��riv��es partielles elliptique aux coefficients variables portant sur la pression. Ont ��t�� ainsi d��velopp��es une m��thode multigrille g��om��trique employant une m��thode de construction de grille grossi��re alg��brique et une variante �� base de semi-grossissement flexible. La strat��gie d'acc��l��ration de convergence de la m��thode fortement coupl��e s'appuie sur deux pivots essentiels : une acc��l��ration interne consistant �� r��soudre chaque syst��me lin��aire par m��thode de Krylov munie de pr��conditionnement par blocs, une acc��l��ration externe par un recours soit �� un sch��ma multigrille non-lin��aire soit �� une m��thode �� sous-espace. Les r��sultats obtenus d��montrent la robustesse des techniques d'acc��l��ration de convergence propos��es. A ��galement ��t�� prouv��e la capacit�� de l'approche fortement coupl��e �� pr��dire de fa��on pr��cise des ��coulements laminaires ou turbulents autour de g��om��tries mod��r��ment complexes en un temps de calcul minime., Th��se de doctorat soutenue le 16/11/1998. @phdthesis{vass:98, title = {Etude num\'erique de techniques d'acc\'el\'eration de convergence lors de la r\'esolution des \'equations de Navier-Stokes en formulation d\'ecoupl\'ee ou fortement coupl\'ee}, author = {Xavier Vasseur}, school = {Universit\'e de Nantes, Ecole Centrale de Nantes, France}, year = {1998} }

In this work, we propose another alternative to find an initial velocity model for the acoustic FWI without any physical knowledge. Motivated by the recent growth of high performance computing (HPC), we tackle the high non-linearity of the problem to minimize, using global optimization methods which are easy to parallelize, in particular, evolution strategies. The first contribution adapt evolution strategies to the FWI setting where the cost function evaluation is the most expensive part. The second contribution is the parameterization of the regarded problem, by being able to represent the model, as faithfully as possible, while limiting the number of parameters needed, since each additional parameter is an additional dimension to explore. The last contribution is to propose a highly parallel evolution strategy adapted to the FWI setting. The initial results on the Salt Dome velocity model using low frequency range, show that great improvement can be done to automate the FWI.