Second order adjoints for solving PDE-constrained optimization problems

Inverse problems are of utmost importance in many fields of science and engineering. In the variational approach inverse problems are formulated as PDE-constrained optimization problems, where the optimal estimate of the uncertain parameters is the minimizer of a certain cost functional subject to the constraints posed by the model equations. The numerical solution of such optimization problems requires the computation of derivatives of the model output with respect to model parameters. The first order derivatives of a cost functional (defined on the model output) with respect to a large number of model parameters can be calculated efficiently through first order adjoint sensitivity analysis. Second order adjoint models give second derivative information in the form of matrix-vector products between the Hessian of the cost functional and user defined vectors. Traditionally, the construction of second order derivatives for large scale models has been considered too costly. Consequently, data assimilation applications employ optimization algorithms that use only first order derivative information, like nonlinear conjugate gradients and quasi-Newton methods. In this paper we discuss the mathematical foundations of second order adjoint sensitivity analysis and show that it provides an efficient approach to obtain Hessian-vector products. We study the benefits of using of second order information in the numerical optimization process for data assimilation applications. The numerical studies are performed in a twin experiment setting with a two-dimensional shallow water model. Different scenarios are considered with different discretization approaches, observation sets, and noise levels. Optimization algorithms that employ second order derivatives are tested against widely used methods that require only first order derivatives. Conclusions are drawn regarding the potential benefits and the limitations of using high-order information in large scale data assimilation problems

Factorization in block-triangularly implicit methods for shallow water applications

The system of first-order ordinary differential equations obtained by spatial discretization of the initial-boundary value problems modelling phenomena in shallow water in 3 spatial dimensions have righthand sides of the form f(t,y) := f1(t,y) + f2(t,y) + f3(t,y) + f4(t,y), where f1, f2 and f3 contain the spatial derivative terms with respect to the x, y and z directions, respectively, and f4 represents the forcing terms and/or reaction terms. The number N of components of f is usually extremely large. It is typical for shallow water applications that the function f4 is nonstiff and that the function f3 corresponding with the vertical spatial direction is much more stiff than the functions f1 and f2 corresponding with the horizontal spatial directions. The reason is that in shallow seas the gridsize in the vertical direction is several orders of magnitude smaller than in the horizontal directions. In order to solve the initial value problem for the system of ordinary differential equations numerically, we need a stiff solver. Stiff IVP solvers are necessarily implicit, requiring the solution of large systems of implicit relations. In a few earlier papers, we considered implicit Runge-Kutta methods leading to fully coupled, implicit systems whose dimension is a multiple of N, and block-diagonally implicit methods in which the implicit relations can be decoupled into subsystems of dimension N. In the present paper, we analyse Rosenbrock type methods and the related DIRK methods (diagonally implicit Runge-Kutta methods) leading to block-triangularly implicit relations. In particular, we shall present a convergence analysis of various iterative methods based on approximate factorization for solving the triangularly implicit relations

Approximate factorization for time-dependent partial differential equations

AbstractThe first application of approximate factorization in the numerical solution of time-dependent partial differential equations (PDEs) can be traced back to the celebrated papers of Peaceman and Rachford and of Douglas of 1955. For linear problems, the Peaceman–Rachford–Douglas method can be derived from the Crank–Nicolson method by the approximate factorization of the system matrix in the linear system to be solved. This factorization is based on a splitting of the system matrix. In the numerical solution of time-dependent PDEs we often encounter linear systems whose system matrix has a complicated structure, but can be split into a sum of matrices with a simple structure. In such cases, it is attractive to replace the system matrix by an approximate factorization based on this splitting. This contribution surveys various possibilities for applying approximate factorization to PDEs and presents a number of new stability results for the resulting integration methods

A comparison of operator splitting and approximate matrix factorization for the shallow water equations in spherical geometry

spherical geometry;The shallow water equations (SWEs) in spherical geometry provide abasic prototypefor developing and testing numerical algorithms for solving the horizontaldynamics in global atmospheric circulation models. When solving the SWEs on a global fine uniform lat-lon grid, an explicit timeintegration method suffers from a severe stability restriction on theadmissible step size. In a previous paper, we investigated an A-stable,linearly-implicit, third-order time integration method (Ros3), which wecombinedwith approximate matrix factorization (AMF) to make it cost-effective. Inthis paper, we further explore this method and we compare itto a Strang-type operator splitting method. Our main focus is on the localerror of the methods, their numerical dispersion relation and their accuracyand efficiency when applied to the well-known SWEs test set. Thecomparison shows that Ros3with AMF accurately presents both low and mid frequency waves. Moreover,Ros3 with AMF makes a good candidate for theefficient solution of the SWEs on a global fine lat-longrid. In contrast, Strang splitting is not advocated, in view of itsinaccuracy in the polar area and the resulting inefficiency

Discontinuous Galerkin methods for solving the acoustic wave equation

In this work we develop a numerical simulator for the propagation of elastic waves by solving the one-dimensional acoustic wave equation with Absorbing Boundary Conditions (ABC’s) on the computational boundaries using Discontinuous Galerkin Finite Element Methods (DGFEM). The DGFEM allows us to easily simulate the presence of a fracture in the elastic medium by means of a linear-slip model. We analize the behaviour of our algorithm by comparing its results against analytic solutions. Furthermore, we show the frequency-dependent effect on the propagation produced by the fracture as appears in previous works. Finally, we present an analysis of the numerical parameters of the method.Fil: Castromán, Gabriel Alejandro. Consejo Nacional de Investigaciones Científicas y Técnicas; Argentina. Universidad Nacional de la Plata. Facultad de Ciencias Astronómicas y Geofísicas. Departamento de Geofísica Aplicada; ArgentinaFil: Zyserman, Fabio Ivan. Consejo Nacional de Investigaciones Científicas y Técnicas; Argentina. Universidad Nacional de la Plata. Facultad de Ciencias Astronómicas y Geofísicas. Departamento de Geofísica Aplicada; Argentin