On the evaluation of general sparse hybrid linear solvers

General sparse hybrid solvers are commonly used kernels for solving wide range of scientific and engineering problems. This work addresses the current problems of efficiently solving general sparse linear equations with direct/iterative hybrid solvers on many core distributed clusters. We briefly discuss the solution stages of Maphys, HIPS, and PDSLin hybrid solvers for large sparse linear systems with their major algorithmic differences. In this category of solvers, different methods with sophisticated preconditioning algorithms are suggested to solve the trade off between memory and convergence. Such solutions require a certain hierarchical level of parallelism more suitable for modern supercomputers that allow to scale for thousand numbers of processors using Schur complement framework. We study the effect of reordering and analyze the performance, scalability as well as memory for each solve phase of PDSLin, Maphys, and HIPS hybrid solvers using large set of challenging matrices arising from different actual applications and compare the results with SuperLU_DIST direct solver. We specifically focus on the level of parallel mechanisms used by the hybrid solvers and the effect on scalability. Tuning and Analysis Utilities (TAU) is employed to assess the efficient usage of heap memory profile and measuring communication volume. The tests are run on high performance large memory clusters using up to 512 processors.     

A block preconditioning technique for the streamfunction‐vorticity formulation of the Navier‐Stokes equations

Iterative methods of Krylov‐subspace type can be very effective solvers for matrix systems resulting from partial differential equations if appropriate preconditioning is employed. We describe and test block preconditioners based on a Schur complement approximation which uses a multigrid method for finite element approximations of the linearized incompressible Navier‐Stokes equations in streamfunction and vorticity formulation. By using a Picard iteration, we use this technology to solve fully nonlinear Navier‐Stokes problems. The solvers which result scale very well with problem parameters. © 2011 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2011     

Dual-primal isogeometric tearing and interconnecting solvers for multipatch continuous and discontinuous Galerkin IgA equations

We present results regarding fast and robust solvers for equations arising from continuous and discontinuous Galerkin discretization of heterogeneous diffusion problems in the context of Isogeometric Analysis. The solvers considered belong to the class of non-overlapping domain decomposition methods which use tearing and interconnection strategies. (© 2016 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Fast Fourier Transform Solvers and Preconditioners for Quadratic Spline Collocation

Quadratic Spline Collocation (QSC) methods of optimal order of convergence have been recently developed for the solution of elliptic Partial Differential Equations (PDEs). In this paper, linear solvers based on Fast Fourier Transforms (FFT)are developed for the solution of the QSC equations. The complexity of the FFT solvers is O(N ² log N), where N is the gridsize in one dimension. These direct solvers can handle PDEs with coefficients in one variable or constant, and Dirichlet, Neumann, alternating Dirichlet-Neumann or periodic boundary conditions, along at least one direction of a rectangular domain. General variable coefficient PDEs are handled by preconditioned iterative solvers. The preconditioner is the QSC matrix arising from a constant coefficient PDE. The convergence analysis of the preconditioner is presented. It is shown that, under certain conditions, the convergence rate is independent of the gridsize. The preconditioner is solved by FFT techniques, and integrated with one-step or acceleration methods, giving rise to asymptotically almost optimal linear solvers, with complexity O(N ² log N). Numerical experiments verify the effectiveness of the solvers and preconditioners, even on problems more general than the analysis assumes. The development and analysis of FFT solvers and preconditioners is extended to QSC equations corresponding to systems of elliptic PDEs.     

Compress‐and‐restart block Krylov subspace methods for Sylvester matrix equations

Block Krylov subspace methods (KSMs) comprise building blocks in many state‐of‐the‐art solvers for large‐scale matrix equations as they arise, for example, from the discretization of partial differential equations. While extended and rational block Krylov subspace methods provide a major reduction in iteration counts over polynomial block KSMs, they also require reliable solvers for the coefficient matrices, and these solvers are often iterative methods themselves. It is not hard to devise scenarios in which the available memory, and consequently the dimension of the Krylov subspace, is limited. In such scenarios for linear systems and eigenvalue problems, restarting is a well‐explored technique for mitigating memory constraints. In this work, such restarting techniques are applied to polynomial KSMs for matrix equations with a compression step to control the growing rank of the residual. An error analysis is also performed, leading to heuristics for dynamically adjusting the basis size in each restart cycle. A panel of numerical experiments demonstrates the effectiveness of the new method with respect to extended block KSMs.     

The relaxation schemes for systems of conservation laws in arbitrary space dimensions

We present a class of numerical schemes (called the relaxation schemes) for systems of conservation laws in several space dimensions. The idea is to use a local relaxation approximation. We construct a linear hyperbolic system with a stiff lower order term that approximates the original system with a small dissipative correction. The new system can be solved by underresolved stable numerical discretizations without using either Riemann solvers spatially or a nonlinear system of algebraic equations solvers temporally. Numerical results for 1-D and 2-D problems are presented. The second-order schemes are shown to be total variation diminishing (TVD) in the zero relaxation limit for scalar equations. ©1995 John Wiley & Sons, Inc.     

Robust Schur complement method for strongly anisotropic elliptic equations

We present robust and asymptotically optimal iterative methods for solving 2D anisotropic elliptic equations with strongly jumping coefficients, where the direction of anisotropy may change sharply between adjacent subdomains. The idea of a stable preconditioning for the Schur complement matrix is based on the use of an exotic non‐conformal coarse mesh space and on a special clustering of the edge space components according to the anisotropy behavior. Our method extends the well known BPS interface preconditioner [2] to the case of anisotropic equations. The technique proposed also provides robust solvers for isotropic equations in the presence of degenerate geometries, in particular, in domains composed of thin substructures. Numerical experiments confirm efficiency and robustness of the algorithms for the complicated problems with strongly varying diffusion and anisotropy coefficients as well as for the isotropic diffusion equations in the ‘brick and mortar’ structures involving subdomains with high aspect ratios. Copyright © 1999 John Wiley & Sons, Ltd.     

Integral equations requiring small numbers of Krylov-subspace iterations for two-dimensional smooth penetrable scattering problems

This paper presents a class of boundary integral equations for the solution of problems of electromagnetic and acoustic scattering by two-dimensional homogeneous penetrable scatterers with smooth boundaries. The new integral equations, which, as is established in this paper, are uniquely solvable Fredholm equations of the second kind, result from representations of fields as combinations of single and double layer potentials acting on appropriately chosen regularizing operators. As demonstrated in this text by means of a variety of numerical examples (that resulted from a high-order Nyström computational implementation of the new equations), these “regularized combined equations” can give rise to important reductions in computational costs, for a given accuracy, over those resulting from previous iterative boundary integral equation solvers for transmission problems.     

Coupling methods for heat transfer and heat flow: Operator splitting and the parareal algorithm

We propose an operator splitting method for coupling heat transfer and heat flow equations. This work is motivated by the need to couple independent industrial heat transfer solvers (e.g., the Aura-Fluid software package) and heat flow solvers (e.g., Openfoam). Such packages are often used to simulate the influence of solar heat in car bodies and are coupled by A-B splitting techniques. The main goal of this work is the acceleration of the coupled software system by iterative operator splitting methods and additional time-parallelism using the parareal algorithm. We present these new splitting techniques along with some novel convergence results and test the splitting-parareal combination on various numerical problems.     

Coupled Models and Parallel Simulations for Three-Dimensional Full-Stokes Ice Sheet Modeling

A three-dimensional full-Stokes computational model is considered for determining the dynamics,temperature,and thickness of ice sheets.The goveming thermomechanical equations consist of the three-dimensional full-Stokes system with nonlinear rheology for the momentum,an advective-diffusion energy equation for temperature evolution,and a mass conservation equation for ice-thickness changes.Here,we discuss the variable resolution meshes,the finite element discretizations,and the parallel algorithms employed by the model components.The solvers are integrated through a well-designed coupler for the exchange of parametric data between components.The discretization utilizes high-quality,variable-resolution centroidal Voronoi Delaunay triangulation meshing and existing parallel solvers.We demonstrate the gridding technology,discretization schemes,and the efficiency and scalability of the parallel solvers through computational experiments using both simplified geometries arising from benchmark test problems and a realistic Greenland ice sheet geometry.     

Spectral discretizations of 3-D elliptic problems and fast domain decomposition methods

An important for applications, the class of hp discretizations of second-order elliptic equations consists of discretizations based on spectral finite elements. The development of fast domain decomposition algorithms for them was restrained by the absence of fast solvers for the basic components of the method, i.e., for local interior problems on decomposition subdomains and their faces. Recently, the authors have established that such solvers can be designed using special factorized preconditioners. In turn, factorized preconditioners are constructed using an important analogy between the stiffness matrices of spectral and hierarchical basis hp-elements (coordinate functions of the latter are defined as tensor products of integrated Legendre polynomials). Due to this analogy, for matrices of spectral elements, fast solvers can be developed that are similar to those for matrices of hierarchical elements. Based on these facts and previous results on the preconditioning of other components, fast domain decomposition algorithms for spectral discretizations are obtained.     

<Emphasis FontCategory="NonProportional">pyomo.dae</Emphasis>: a modeling and automatic discretization framework for optimization with differential and algebraic equations

We describe pyomo.dae, an open source Python-based modeling framework that enables high-level abstract specification of optimization problems with differential and algebraic equations. The pyomo.dae framework is integrated with the Pyomo open source algebraic modeling language, and is available at http://www.pyomo.org. One key feature of pyomo.dae is that it does not restrict users to standard, predefined forms of differential equations, providing a high degree of modeling flexibility and the ability to express constraints that cannot be easily specified in other modeling frameworks. Other key features of pyomo.dae are the ability to specify optimization problems with high-order differential equations and partial differential equations, defined on restricted domain types, and the ability to automatically transform high-level abstract models into finite-dimensional algebraic problems that can be solved with off-the-shelf solvers. Moreover, pyomo.dae users can leverage existing capabilities of Pyomo to embed differential equation models within stochastic and integer programming models and mathematical programs with equilibrium constraint formulations. Collectively, these features enable the exploration of new modeling concepts, discretization schemes, and the benchmarking of state-of-the-art optimization solvers.     

MIP-based instantaneous control of mixed-integer PDE-constrained gas transport problems

We study the transient optimization of gas transport networks including both discrete controls due to switching of controllable elements and nonlinear fluid dynamics described by the system of isothermal Euler equations, which are partial differential equations in time and 1-dimensional space. This combination leads to mixed-integer optimization problems subject to nonlinear hyperbolic partial differential equations on a graph. We propose an instantaneous control approach in which suitable Euler discretizations yield systems of ordinary differential equations on a graph. This networked system of ordinary differential equations is shown to be well-posed and affine-linear solutions of these systems are derived analytically. As a consequence, finite-dimensional mixed-integer linear optimization problems are obtained for every time step that can be solved to global optimality using general-purpose solvers. We illustrate our approach in practice by presenting numerical results on a realistic gas transport network.     

Parallel preconditioners and multigrid solvers for stochastic polynomial chaos discretizations of the diffusion equation at the large scale

This paper presents parallel preconditioners and multigrid solvers for solving linear systems of equations arising from stochastic polynomial chaos formulations of the diffusion equation with random coefficients. These preconditioners and solvers are extensions of the preconditioner developed in an earlier paper for strongly coupled systems of elliptic partial differential equations that are norm equivalent to systems that can be factored into an algebraic coupling component and a diagonal differential component. The first preconditioner, which is applied to the norm equivalent system, is obtained by sparsifying the inverse of the algebraic coupling component. This sparsification leads to an efficient method for solving these systems at the large scale, even for problems with large statistical variations in the random coefficients. An extension of this preconditioner leads to stand‐alone multigrid methods that can be applied directly to the actual system rather than to the norm equivalent system. These multigrid methods exploit the algebraic/differential factorization of the norm equivalent systems to produce variable‐decoupled systems on the coarse levels. Moreover, the structure of these methods allows easy software implementation through re‐use of robust high‐performance software such as the Hypre library package. Two‐grid matrix bounds will be established, and numerical results will be given. Copyright © 2015 John Wiley & Sons, Ltd.     

Solving Differential Matrix Equations using Parareal

Differential matrix equations appear in many applications like optimal control of partial differential equations, balanced truncation model order reduction of linear time varying systems and many more. Here, we will focus on differential Riccati equations (DRE). Solving such matrix-valued ordinary differential equations (ODE) is a highly time consuming process. We present a Parareal based algorithm applied to Rosenbrock methods for the solution of the matrix-valued differential Riccati equations. Considering problems of moderate size, direct matrix equation solvers for the solution of the algebraic Lyapunov equations arising inside the time intgration methods are used. (© 2016 Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim)     

A Test Set for stiff Initial Value Problem Solvers in the open source software R: Package deTestSet

In this paper we present the R package deTestSet that includes challenging test problems written as ordinary differential equations (ODEs), differential algebraic equations (DAEs) of index up to 3 and implicit differential equations (IDEs). In addition it includes 6 new codes to solve initial value problems (IVPs). The R package is derived from the Test Set for Initial Value Problem Solvers available at http://www.dm.uniba.it/~testset which includes documentation of the test problems, experimental results from a number of proven solvers, and Fortran subroutines providing a common interface to the defining problem functions. Many of these facilities are now available in the R package deTestSet, which comprises an R interface to the test problems and to most of the Fortran solvers. The package deTestSet is free software which is distributed under the GNU General Public License, as part of the R open source software project.     

Variable-step variable-order 3-stage Hermite-Birkhoff-Obrechkoff DDE solver of order 4 to 14

This article presents a solver for delay differential equations (DDEs) called HBO414DDE based on a hybrid variable-step variable-order 3-stage Hermite-Birkhoff-Obrechkoff ODE solver of order 4 to 14. The current version of our method solves DDEs with state dependent, non-vanishing, small, vanishing and asymptotically vanishing delays, except neutral type and initial value DDEs. Delayed values are computed using Hermite interpolation, small delays are dealt with by extrapolation, and discontinuities are located by a bisection method. HBO414DDE was tested on several problems and results were compared with those of known solvers like SYSDEL and the recent Matlab DDE solver ddesd and statistics show that it gives, most of the time, a smaller relative error than the other solvers for the same number of function evaluations.     

Waveform relaxation methods for implicit differential equations

We apply a Runge-Kutta-based waveform relaxation method to initial-value problems for implicit differential equations. In the implementation of such methods, a sequence of nonlinear systems has to be solved iteratively in each step of the integration process. The size of these systems increases linearly with the number of stages of the underlying Runge-Kutta method, resulting in high linear algebra costs in the iterative process for high-order Runge-Kutta methods. In our earlier investigations of iterative solvers for implicit initial-value problems, we designed an iteration method in which the linear algebra costs are almost independent of the number of stages when implemented on a parallel computer system. In this paper, we use this parallel iteration process in the Runge-Kutta waveform relaxation method. In particular, we analyse the convergence of the method. The theoretical results are illustrated by a few numerical examples.     

A Lagrangian Scheme à la Brenier for the Incompressible Euler Equations

We approximate the regular solutions of the incompressible Euler equations by the solution of ODEs on finite-dimensional spaces. Our approach combines Arnold’s interpretation of the solution of the Euler equations for incompressible and inviscid fluids as geodesics in the space of measure-preserving diffeomorphisms, and an extrinsic approximation of the equations of geodesics due to Brenier. Using recently developed semi-discrete optimal transport solvers, this approach yields a numerical scheme which is able to handle problems of realistic size in 2D. Our purpose in this article is to establish the convergence of this scheme towards regular solutions of the incompressible Euler equations, and to provide numerical experiments on a few simple test cases in 2D.     

Quasi-Newton acceleration for equality-constrained minimization

Optimality (or KKT) systems arise as primal-dual stationarity conditions for constrained optimization problems. Under suitable constraint qualifications, local minimizers satisfy KKT equations but, unfortunately, many other stationary points (including, perhaps, maximizers) may solve these nonlinear systems too. For this reason, nonlinear-programming solvers make strong use of the minimization structure and the naive use of nonlinear-system solvers in optimization may lead to spurious solutions. Nevertheless, in the basin of attraction of a minimizer, nonlinear-system solvers may be quite efficient. In this paper quasi-Newton methods for solving nonlinear systems are used as accelerators of nonlinear-programming (augmented Lagrangian) algorithms, with equality constraints. A periodically-restarted memoryless symmetric rank-one (SR1) correction method is introduced for that purpose. Convergence results are given and numerical experiments that confirm that the acceleration is effective are presented. This work was supported by FAPESP, CNPq, PRONEX-Optimization (CNPq / FAPERJ), FAEPEX, UNICAMP.