A massively parallel exponential integrator for advection-diffusion models

This work considers the Real Leja Points Method (ReLPM) for the exponential integration of large-scale sparse systems of ODEs, generated by Finite Element or Finite Difference discretizations of 3-D advection-diffusion models. We present an efficient parallel implementation of ReLPM for polynomial interpolation of the matrix exponential propagators. A scalability analysis of the most important computational kernel inside the code, the parallel sparse matrix\u2013vector product, has been performed, as well as an experimental study of the communication overhead. As a result of this study an optimized parallel sparse matrix\u2013vector product routine has been implemented. The resulting code shows good scaling behavior even when using more than one thousand processors. The numerical results presented on a number of very large test cases gives experimental evidence that ReLPM is a reliable and efficient tool for the simulation of complex hydrodynamic processes on parallel architectures

Optimal control and robust estimation for ocean wave energy converters

This thesis deals with the optimal control of wave energy converters and some associated observer design problems. The first part of the thesis will investigate model predictive control of an ocean wave energy converter to maximize extracted power. A generic heaving converter that can have both linear dampers and active elements as a power take-off system is considered and an efficient optimal control algorithm is developed for use within a receding horizon control framework. The optimal control is also characterized analytically. A direct transcription of the optimal control problem is also considered as a general nonlinear program. A variation of the projected gradient optimization scheme is formulated and shown to be feasible and computationally inexpensive compared to a standard nonlinear program solver. Since the system model is bilinear and the cost function is not convex quadratic, the resulting optimization problem is shown not to be a quadratic program. Results are compared with other methods like optimal latching to demonstrate the improvement in absorbed power under irregular sea condition simulations. In the second part, robust estimation of the radiation forces and states inherent in the optimal control of wave energy converters is considered. Motivated by this, low order H∞ observer design for bilinear systems with input constraints is investigated and numerically tractable methods for design are developed. A bilinear Luenberger type observer is formulated and the resulting synthesis problem reformulated as that for a linear parameter varying system. A bilinear matrix inequality problem is then solved to find nominal and robust quadratically stable observers. The performance of these observers is compared with that of an extended Kalman filter. The robustness of the observers to parameter uncertainty and to variation in the radiation subsystem model order is also investigated. This thesis also explores the numerical integration of bilinear control systems with zero-order hold on the control inputs. Making use of exponential integrators, exact to high accuracy integration is proposed for such systems. New a priori bounds are derived on the computational complexity of integrating bilinear systems with a given error tolerance. Employing our new bounds on computational complexity, we propose a direct exponential integrator to solve bilinear ODEs via the solution of sparse linear systems of equations. Based on this, a novel sparse direct collocation of bilinear systems for optimal control is proposed. These integration schemes are also used within the indirect optimal control method discussed in the first part.Open Acces

A black-box rational Arnoldi variant for Cauchy-Stieltjes matrix functions

Rational Arnoldi is a powerful method for approximating functions of large sparse matrices times a vector. The selection of asymptotically optimal parameters for this method is crucial for its fast convergence. We present and investigate a novel strategy for the automated parameter selection when the function to be approximated is of Cauchy-Stieltjes (or Markov) type, such as the matrix square root or the logarithm. The performance of this approach is demonstrated by numerical examples involving symmetric and nonsymmetric matrices. These examples suggest that our black-box method performs at least as well, and typically better, as the standard rational Arnoldi method with parameters being manually optimized for a given matrix

Matrices, Moments and Quadrature: Applications to Time- Dependent Partial Differential Equations

The numerical solution of a time-dependent PDE generally involves the solution of a stiff system of ODEs arising from spatial discretization of the PDE. There are many methods in the literature for solving such systems, such as exponential propagation iterative (EPI) methods, that rely on Krylov projection to compute matrix function-vector products. Unfortunately, as spatial resolution increases, these products require an increasing number of Krylov projection steps, thus drastically increasing computational expense

Approximation of the Neutron Diffusion Equation on Hexagonal Geometries

La ecuación de la difusión neutrónica describe la población de neutrones de un reactor nuclear. Este trabajo trata con este modelo para reactores nucleares con geometría hexagonal. En primer lugar se estudia la ecuación de la difusión neutrónica. Este es un problema diferencial de valores propios, llamado problema de los modos Lambda. Para resolver el problema de los modos Lambda se han comparado diferentes métodos en geometrías unidimensionales, resultando como el mejor el método de elementos espectrales. Usando este método discretizamos los operadores en geometrías bidimensiones y tridimensionales, resolviendo el problema algebraica de valores propios resultante con el método de Arnoldi. La distribución de neutrones estado estacionario se utiliza como condición inicial para la integración de la ecuación de la difusión neutrónica dependiente del tiempo. Se utiliza un método de Euler implícito para integrar en el tiempo. Cuando un nodo está parcialmente insertado aparece un comportamiento no físico de la solución, el efecto ``rod cusping'', que se corrige mediante la ponderación de las secciones eficaces con el flujo del paso de tiempo anterior. Cuando la solución de los sistemas algebraicos que surgen en el método hacia atrás, un método de Krylov se utiliza para resolver los sistemas resultantes, y diferentes estrategias de precondicionamiento se evalúan se. La primera consiste en el uso de la estructura de bloque obtenido por los grupos de energía para resolver el sistema por bloques, y diferentes técnicas de aceleración para el esquema iterativo de bloques y un precondicionador utilizando esta estructura de bloque se proponen. Además se estudia un precondicionador espectral, que hace uso de la información en un subespacio de Krylov para precondicionar el siguiente sistema. También se proponen métodos exponenciales de segundo y cuarto orden integrar la ecuación de difusión neutrónica dependiente del tiempo, donde la exponencial de la matriz del sistema tiene quGonzález Pintor, S. (2012). Approximation of the Neutron Diffusion Equation on Hexagonal Geometries [Tesis doctoral no publicada]. Universitat Politècnica de València. https://doi.org/10.4995/Thesis/10251/17829Palanci

Improving Efficiency of Rational Krylov Subspace Methods

This thesis studies two classes of numerical linear algebra problems, approximating the product of a function of a matrix with a vector, and solving the linear eigenvalue problem $Av=\lambda Bv$ for a small number of eigenvalues. These problems are solved by rational Krylov subspace methods (RKSM). We present several improvements in two directions: pole selection and applying inexact methods. In Chapter 3, a flexible extended Krylov subspace method ($\mathcal{F}$-EKSM) is considered for numerical approximation of the action of a matrix function $f(A)$ to a vector $b$, where the function $f$ is of Markov type. $\mathcal{F}$-EKSM has the same framework as the extended Krylov subspace method (EKSM), replacing the zero pole in EKSM with a properly chosen fixed nonzero poles. For symmetric positive definite matrices, the optimal fixed pole is derived for $\mathcal{F}$-EKSM to achieve the lowest possible upper bound on the asymptotic convergence factor, which is lower than that of EKSM. The analysis is based on properties of Faber polynomials of $A$ and $(I-A/s)^{-1}$. For large and sparse matrices that can be handled efficiently by LU factorizations, numerical experiments show that $\mathcal{F}$-EKSM and a variant of RKSM based on a small number of fixed poles outperform EKSM in both storage and runtime, and they usually have advantage over adaptive RKSM in runtime. Chapter 4 concerns the theory and development of inexact RKSM for approximating the action of a function of matrix $f(A)$ to a column vector $b$. At each step of RKSM, a shifted linear system of equations needs to be solved to enlarge the subspace. For large-scale problems, arising from discretizations of PDEs in 3D domains, such a linear system is usually solved by an iterative method approximately. The main question is how to relax the accuracy of these linear solves without negatively affecting the convergence for approximating $f(A)b$. Our insight into this issue is obtained by exploring the residual bounds on the rational Krylov subspace approximations to $f(A)b$, based on the decaying behavior of the entries in the first column of the matrix function of the block Rayleigh quotient of $A$ with respect to the rational Krylov subspaces. The decay bounds on these entries for both analytic functions and Markov functions can be efficiently and accurately evaluated by appropriate quadrature rules. A heuristic based on these bounds is proposed to relax the tolerances of the linear solves arising from each step of RKSM. As the algorithm progresses toward convergence, the linear solves can be performed with increasingly lower accuracy and computational cost. Numerical experiments for large nonsymmetric matrices show the effectiveness of the tolerance relaxation strategy for the inexact linear solves of RKSM. In Chapter 5, inexact RKSM are studied to solve large-scale nonsymmetric eigenvalue problems. Similar to the problem setting in Chapter 4, each iteration (outer step) of RKSM requires solution to a shifted linear system to enlarge the subspace, but these linear solves by direct methods are prohibitive due to the problem scale. Errors are introduced at each outer step if these linear systems are solved approximately by iterative methods (inner step), and these errors accumulate in the rational Krylov subspace. In this thesis, we derive an upper bound on the errors that can be introduced at each outer step to maintain the same convergence as exact RKSM for approximating an invariant subspace. Since this bound is inversely proportional to the current eigenresidual norm of the desired invariant subspace, the tolerance of iterative linear solves at each outer step can be relaxed with the outer iteration progress. A restarted variant of the inexact RKSM is also proposed. Numerical experiments show the effectiveness of relaxing the inner tolerance to save computational cost

Investigation on polynomial integrators for time-domain electromagnetics using a high-order discontinuous Galerkin method

International audienceIn this work, we investigate the application of polynomial integrators in a high-order discontinuous Galerkin method for solving the time-domain Maxwell equations. After the spatial discretization, we obtain a time-continuous system of ordinary differential equations of the form, ∂tY(t)=HY(t), where Y(t) is the electromagnetic field, H is a matrix containing the spatial derivatives, and t is the time variable. The formal solution is written as the exponential evolution operator, exp(tH), acting on a vector representing the initial condition of the electromagnetic field. The polynomial integrators are based on the approximation of exp(tH) by an expansion of the form ∑ _m=0^\infinity gm(t) Pm(H), where gm(t) is a function of time and Pm(H) is a polynomial of order m satisfying a short recursion. We introduce a general family of expansions of exp(tH) based on Faber polynomials. This family of expansions is suitable for non-Hermitian matrices, and consequently the proposed integrators can handle absorbing media and conductive materials. We discuss the efficient implementation of this technique, and based on some test problems, we compare the virtues and shortcomings of the algorithm. We also demonstrate how this scheme provides an efficient alternative to standard explicit integrators