CM2018: FIFTEENTH COPPER MOUNTAIN CONFERENCE ON ITERATIVE METHODS
PROGRAM FOR FRIDAY, MARCH 30TH
Days:
previous day
all days

View: session overviewtalk overview

08:00-10:05 Session 14A: Multilevel Solvers. Room: Bighorn B.
08:00
Two-level Fourier analysis of multigrid for higher-order finite-element methods

ABSTRACT. Local Fourier analysis (LFA) is well-known as a predictive design and analysis tool for multigrid methods and has been applied to many problems in many contexts. Here, we consider the application of LFA to analyze the convergence properties of multigrid methods for higher-order finite-element approximations to the Laplacian problem. We find that the LFA smoothing factor fails to predict the true performance of multigrid in this setting. This failure is explained, and we propose an alternative analysis that yields a reasonable prediction of performance that can be used, for example, to help choose the correct damping parameters for relaxation. Finally, we present two-grid and multigrid experiments, demonstrating that the corrected parameter choice yields a significant improvement in the resulting convergence factors.

08:25
Improving the Performance of Line Relaxation in BoxMG

ABSTRACT. Line and plane relaxation are important components for robust variational multigrid methods on structured grids when using standard coarsening with anisotropic problems. The key computational element for these relaxation methods is the solution of tridiagonal systems of linear equations. As line and plane relaxation are typically the dominant cost for these multigrid methods, an efficient tridiagonal solver is key to their performance. The commonly used Thomas algorithm is strictly serial so alternative algorithms such as Wang's two-level algorithm, and cyclic reduction have been explored. In this talk, we present the results of extending Wang's algorithm to a multilevel setting to improve the performance of the tridiagonal solver in Black Box Multigrid (BoxMG).

08:50
LFA Lab - Local Fourier Analysis including periodic stencils in Python

ABSTRACT. Local Fourier Analysis (LFA) is an established tool for the analysis of multigrid methods and similar iterative methods. Recently, we introduced a generalization for periodic stencils, allowing for the analysis of non-constant coefficient problems and block smoothers. In order to ease the analysis using our methodology we developed a new LFA tool usable in Python partly implemented in C++.

In this talk LFA including our extension will be briefly reviewed and our tool will be presented.

09:15
A Posteriori Error Estimators for the Frank-Oseen Liquid Crystal Model
SPEAKER: David Emerson

ABSTRACT. A posteriori error estimators aim to provide easily computable and reliable bounds on the error of numerical solutions for partial differential equations (PDEs) and variational systems. Accurate error estimators significantly increase the efficiency of numerical methods by facilitating the construction of optimal discretizations via adaptive mesh refinement (AMR). Furthermore, such estimators offer a means of objectively measuring the quality of a computed numerical solution. In this talk, we discuss the construction of novel and reliable a posteriori error estimators for the nonlinear variational system arising in the context of the Frank-Oseen free-energy model of cholesteric liquid crystals. Such estimators are based on a general framework established by Verf\"{u}rth for nonlinear PDEs. The proposed error estimators are composed of readily computable quantities on each element of a finite-element mesh, enabling the formulation of efficient AMR strategies. A set of numerical examples is presented wherein the AMR schemes successfully provide significant reductions in computational work while producing solutions that are highly competitive with those found on uniformly refined meshes. The numerical simulations include error analysis for problems with known analytical solutions as well as performance metrics for configurations for which only numerical solutions are currently available.

09:40
Iterative Algebraic Smoothing Aggregation Multiscale Solver
SPEAKER: Pavel Tomin

ABSTRACT. In recent years, the development of the physics-based preconditioners has received significant attention in the numerical simulation community, to name a few here as fixed stress point method, deflation methods, multiscale methods, etc. Usually, preconditioners are used to damp slowly varying error modes in the linear solver stage of the problem residuals. State-of-the-art multiscale solvers use a sequence of aggressive restriction, coarse-grid correction and prolongation operators to handle low-frequency modes on the coarse grid. High-frequency errors are then resolved by employing a smoother on the original fine grid.

In this paper, an iterative algebraic smoothing aggregation multiscale (iASAM) solver is proposed. The prolongation and restriction operators in this method are constructed by using variables aggregation and Jacobi smoothing iterations of indicator vectors corresponding to these aggregations. The proposed method generalizes the existing multiscale restriction smoothed basis (MsRSB) method.

The new method does not require any coarse partitioning and, hence, can be applied to a general unstructured topology of the fine scale. The iASAM provides a reasonably good approximation and speed up the convergence when used as a preconditioner in an iterative fine-scale solver. In addition, the method obeys a good theoretical scalability essential for parallel simulations. Numerical results for selected test cases (e.g., coupled fluid flow and geomechanics, including geomechanics contact problems) are presented, discussed and future studies are outlined.

08:00-10:05 Session 14B: Multi-physics EM/Fluid Plasma. Room: Bighorn C/1.
08:00
Conjugate gradient for nonsingular saddle-point systems with a highly singular leading block
SPEAKER: unknown

ABSTRACT. We consider iterative solvers for large, sparse, symmetric linear systems with a saddle-point structure. Since such systems are indefinite, the conjugate gradient (CG) method cannot typically be used. However, in the case of a maximally rank-deficient leading block, we show that there are two necessary and sufficient conditions that allow for CG to be used. We show that the conditions are satisfied for a model mixed Maxwell problem. To support our analysis, we present several three-dimensional numerical experiments on complicated computational domains or with variable coefficients.

08:25
A 6th order finite volume multigrid Poisson solver

ABSTRACT. The numerical simulation of electrostatic plasma interactions using a multifluid model [1] requires the coupling of a high-order finite-volume Poisson solver with the fluid solver. In this talk, we present a 6th order finite volume scheme for the numerical solution of Poisson's equation. The starting point of our discussion is the derivation of a novel compact finite volume discretization of the Laplacian which involves a 5^3 stencil. This is a generalization of the L_{7} and L_{19} Laplacian discretizations that have appeared in [2] and [3] and that are 2nd and 4th order accurate, respectively. A suitable preconditioning (Mehrstellen correction) of the right hand side of Poisson's equation, that is essential for obtaining 6th order convergence, is also presented. Following [2] we describe a multigrid extension of our 6th order scheme on a hierarchy of locally refined block structured grids. Handling the computation of the discrete Laplacian at the interfaces between grids at different adaptation levels is accomplished with the interpolation schemes presented in [4]. Several numerical examples are discussed to demonstrate the convergence properties of the method.

[1] D. Ghosh, C. Kavouklis, T. Chapman, R.L. Berger, Numerical Simulation of Counterstreaming Plasma Interactions using a Multifluid Model, 13th World Congress on Computational Mechanics, New York, NY, 2018.

[2] D.F. Martin, K.L. Cartwright, Solving Poisson's Equation using Adaptive Mesh Refinement, EECS Department, Technical Report UCB/ERL, M96/66, University of California, Berkeley, 1996.

[3] M. Barad, P. Colella, A Fourth-Order Accurate Local Refinement Method for Poisson's Equation, Journal of Computational Physics, 2005.

[4] P. McCorquodale, P. Colella, A High-Order Finite-Volume Method for Conservation Laws on Locally Refined Grids, Communications in Applied Mathematics and Computational Science, 2011.

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344 and funded by the LDRD Program at LLNL under project tracking code 17-ERD-081.

08:50
Robust and Scalable Preconditioning of a Newton-Krylov Solver for All-Speed Melt Pool Flow Physics
SPEAKER: unknown

ABSTRACT. In this paper, we introduce a multigrid block-based preconditioner for solving linear systems arising from a Discontinuous Galerkin discretization of the all-speed Navier-Stokes equations with phase change. The equations are discretized in conservative form with a reconstructed Discontinuous Galerkin (rDG) method and integrated with fully-implicit time discretization schemes. To robustly converge the numerically stiff systems, we use the Newton-Krylov framework with a primitive-variable formulation (pressure, velocity, and temperature), which is better conditioned than the conservative-variable form at low-Mach numbers. In the limit of large acoustic CFL number and viscous Fourier number, there is a strong coupling between the velocity-pressure system and the linear systems become non-diagonally dominant, respectively. To effectively solve these ill-conditioned systems, an approximate block factorization preconditioner is developed, which uses the Schur complement to reduce a 3x3 block system into a sequence of two 2x2 block systems: velocity-pressure (vP) and velocity-temperature (vT). We compare the performance of the (vP-vT) Schur complement preconditioner to classic preconditioning strategies: monolithic algebraic multigrid (AMG), element-block SOR, and physics-block Gauss-Seidel. The preconditioned solver performance is investigated in the limit of large CFL and Fourier numbers for compressible lid-driven cavity flow and melt convection of steel. Numerical results demonstrate that the (vP-vT) Schur complement preconditioned solver has excellent algorithmic scalability, parallel scalability, and is robust for highly ill-conditioned systems, for all tested rDG discretization schemes (up to 4th-order).

09:15
A Domain Decomposition Method for efficient and accurate simulations of wave scattering

ABSTRACT. Scattering problems by electrically large objects are of wide interest in many fields. The numerical computation of such problems remains limited by computer resources, due to the large number of unknowns, especially when inhomogeneous materials are present. Domain decomposition methods (DDM) are particularly attractive for the solution of large finite elements problem. It is decomposed into several coupled sub-problems which can be solved independently with accurate transmission conditions, thus reducing considerably the memory storage requirements. Moreover, DDM are intrinsically well suited for numerical implementation on parallel computers. We present here a DDM to solve accurately and efficiently a 3D Maxwell problem, using a coupling between finite elements in the volume and an integral representation (which is an exact radiation condition) on the external boundary $\Gamma_e$. The volume $\Omega_i$ is decomposed into several concentric subdomains $\Omega_i$ (i = 1, ..., N ) and different types of transmission conditions are prescribed on their interfaces. The transmission conditions are then optimized to reduce the condition number of the problem. The global problem (finite elements + integral representation) is solved using a preconditioned GMRES, where the preconditioner is a block diagonal, each block representing either a subdomain $\Omega_i$ or the boundary $\Gamma_e$ . To enhance the efficiency of the solver, a compression algorithm is applied to the resolution of the integral representation, reducing drastically the memory requirement as well as the time of each matrix-vector product. Parallelization is achieved using MPI and threads, for both the solution of the subdomains and the integral representation, and we will present results obtained on the CEA supercomputer Tera 1000.

09:40
On Scalable Block Preconditioning of Multifluid Electromagnetic Plasma Systems*
SPEAKER: John Shadid

ABSTRACT. The mathematical basis for the continuum approximations of multifluid plasma physics systems is the solution of the governing equations describing conservation of mass, momentum, and total energy for each fluid species, along with Maxwell's equations for the electromagnetic fields. The resulting systems are characterized by strong nonlinear and nonsymmetric coupling of fluid and electromagnetic phenomena, as well as the significant range of time- and length-scales that these interactions produce. To enable accurate and stable approximation of these systems a range of spatial and temporal discretization methods are employed. In the context of finite element spatial discretization these can include continuous and discontinuous Galerkin methods of the fluid sub-systems, and divergence cleaning/stabilization or structure-preserving approaches for the electromagnetics sub-system. For effective time integration of the longer time-scale response of these systems some form of implicitness is required. Two well-structured approaches, of recent interest, are fully-implicit and implicit-explicit (IMEX) type methods. The requirement to accommodate disparate spatial discretizations, and allow the flexible assignment of mechanisms as explicit or implicit operators, implies a wide variation in unknown coupling, ordering, the nonzero block structure, and the conditioning of the implicit sub-system. These characteristics make the scalable and efficient iterative solution of these systems extremely challenging.

To enable robust, scalable and efficient solution of the large-scale sparse linear systems generated by a Newton linearization, fully-coupled multilevel preconditioners are developed. The multilevel preconditioners are based on two differing approaches. The first technique employs a graph-based aggregation method applied to the nonzero block structure of the Jacobian matrix [1,2]. The second approach utilizes approximate block factorization (ABF) methods and physics-based preconditioning approaches that reduce the coupled systems into a set of simplified systems to which multilevel methods are applied [3]. To demonstrate the flexibility of implicit/IMEX FE discretizations and the fully-coupled Newton-Krylov-AMG solution approaches various forms of discrete formulations for multifluid electromagnetic plasma models are considered with results presented for representative plasma physics problems of current scientific interest. Additionally, the discussion considers the robustness, efficiency, and the parallel and algorithmic scaling of the preconditioning methods. Weak scaling results include studies on up to 1M cores.

References [1] J. N. Shadid, R. P. Pawlowski, E. C. Cyr, R. S. Tuminaro, L. Chacon, P. D. Weber, Scalable Implicit Incompressible Resistive MHD with Stabilized FE and Fully-coupled Newton-Krylov-AMG, Comput. Methods Appl. Mech. Engrg. 304, 1–25, 2016 [2] P.T. Lin, J.N. Shadid, J.J. Hu, R.P. Pawlowski, E.C. Cyr, "Performance of Fully-coupled Algebraic Multigrid Preconditioners for Large-scale VMS Resistive MHD," Journal of Computational and Applied Math., 2017, in press (https://doi.org/10.1016/j.cam.2017.09.028) [3] E. G. Phillips, J. N. Shadid, E. C. Cyr, Scalable Preconditioners for Structure Preserving Discretizations of Maxwell Equations In First Order Form, submitted to SISC

08:00-10:05 Session 14C: Saddle Point/Indefinite Solvers. Room: Bighorn C/2.
08:00
New stabilized discretization and preconditioner for poroelastic equations
SPEAKER: Peter Ohm

ABSTRACT. In this talk, the popular P1-RT0-P0 discretization of the three-field formulation of Biot’s consolidation problem for porous media will be discussed. In general, this finite-element formulation does not satisfy an inf-sup condition uniformly with respect to the physical parameters and, thus, several issues arise when using it in numerical simulations. Here, we consider a stabilization technique that enriches the piecewise linear finite-element space of the displacement field with the span of edge/face bubble functions. For Biot’s model, this gives rise to discretizations that are uniformly stable with respect to the physical parameters. Additionally, a perturbation of the bilinear form is done, which allows for local elimination of the bubble functions, providing a uniformly stable scheme with the same number of degrees of freedom as the classical P1-RT0-P0 approach. Time permitting, we will also discuss the linear solvers developed for solving the resulting linear systems of such a discretization. More specifically, we will present preconditoners for Krylov iterative methods that are robust with respect to the physical and discretization parameters of the model. Numerical tests confirm the theory for poroelastic test problems.

08:25
A Relaxed Physical Factorization preconditioner for mixed finite element coupled poromechanics.
SPEAKER: unknown

ABSTRACT. In this paper, we introduce a Relaxed Physical Factorization (RPF) preconditioner for the iterative solution of coupled poromechanical problems. The use of mixed finite element method together with implicit time discretization of Biot equations leads to a sequence of unsymmetric and indefinite linear systems with a 3 times 3 block structure. A monolithic approach based on a preconditioned Krylov subspace method is used for efficiently solving such systems. This work focuses on the development and implementation of a new preconditioner inspired by the Relaxed Dimensional Factorization (RDF) preconditioner introduced by Benzi et al. [2011, 2016] for the saddle-point problems arising from Navier-Stokes equations. A similar idea is here used by replacing the Dimensional block matrix subdivision with a Physically-based partitioning. Then, the RDF formulation is extended by taking into account a non-zero contribution from the (3,3) block. A theoretical analysis of the preconditioned matrix eigenspectrum is provided, verifying that it is clustered away from 0. The setup of the proposed preconditioner involves the computation of two inner local preconditioners and the selection of a relaxation parameter 'alpha'. In order to compute the optimal value of 'alpha', a new algorithm is advanced with an automatic control of the possible ill-conditioning of 'alpha'-dependent inner blocks. Two sets of numerical experiments are taken into account. In the former, Mandel's problem, a traditional benchmark in linear poromechanics, is used to check the theoretical properties of the preconditioner. The results validate the accuracy of the 'alpha' estimation independently of the time step size and mesh refinement. In the latter, a real-world consolidation problem is addressed to test the performance and robustness of the proposed algorithm. The results indicate a very fast RPF convergence and a robust behavior versus the nearly optimal relaxation parameter 'alpha'.

M. Benzi, M. Ng, Q. Niu, Z. Wang, A Relaxed Dimensional Factorization preconditioner for the incompressible Navier-Stokes equations, J. Comput. Phys. 230 (16) (2011) 6185-6202.

M. Benzi, S. Deparis, G. Grandperrin, A. Quarteroni, Parameter estimates for the Relaxed Dimensional Factorization preconditioner and application to hemodynamics, Comput. Methods Appl. Mech. Engrg. 300 (2016) 129-145.

08:50
A Tridiagonalization Method for Symmetric and Quasi-Definite Saddle-Point Systems
SPEAKER: Daniel Ruiz

ABSTRACT. We consider the solution of symmetric saddle-point systems with symmetric and positive definite leading block. Such systems arise in numerous applications, including optimization, fluid dynamics, and data assimilation. In the large-scale case, or the case where the blocks are only available as operators, it is common to employ a Krylov method. Prime candidates are the methods MINRES and SYMMLQ of Paige and Saunders, both of which were designed with general symmetric indefinite systems in mind, but neither of which exploits the specific block structure of the system.

The main idea stems from the simple observation that the solution the system with right-hand side (b,c) may be written as the sum of the solutions of the same system with right-hand sides (b,0) and (0,c). In turn, those two systems are the optimality conditions of a least-squares and a least-norm problem.

We propose an approach that meshes an iterative method for least-squares problems with one for least-norm problems in such a way that both problems can be solved in one pass. Each iteration of the proposed procedure has the same cost and almost the same storage requirements as one iteration of MINRES or SYMMLQ but only roughly half as many iterations are required.

Our approach is based on an orthogonal tridiagonalization process initially proposed by Saunders, Simon and Yip in 1988 for square, but not necessarily symmetric, matrices. The process suggests two approaches of the Krylov kind for solving a linear system: one of the minimum-residual kind in the vein of MINRES and named USYMQR, and one of the minimum-norm kind in the vein of SYMMLQ and named USYMLQ. The tridiagonalization process reduces to the symmetric Lanczos process in the symmetric case but differs from the Golub-Kahan biorthogonalization process in that it must be initialized with two vectors. As a consequence, both USYMQR and USYMLQ can be used in conjunction to solve a system and an adjoint system at the same time.

We present preconditioned and regularized variants that apply to symmetric and quasi-definite linear systems and illustrate the performance of our implementation on systems coming from optimization and fluid flow.

We also present variants of the tridiagonalization procedure specifically designed to address rank-deficient problems.

09:15
An inner-outer iterations approach based on the Golub-Kahan bidiagonalization method applied to saddle point problems
SPEAKER: Carola Kruse

ABSTRACT. For the analysis of critical industrial applications, as for example the structural analysis of the containment buildings of nuclear power plants, it is important to have a reliable and efficient simulation software. The open source finite element software Code_Aster developed by EDF (Electricité De France) represents such a tool. So far being only mildly parallel, the goal of the project PAMSIM is to efficiently parallelize the code. Staying with the example of the containment building, we look at a problem in linear elasticity, wherein one-, two- and three-dimensional finite elements (representing the outer shell, metallic wires and the concrete block) are coupled by multi-point constraints (MPCs). The MPCs are enforced by Lagrange multipliers and thus the resulting matrix to solve is indefinite and of a 2x2 block structure, with the (2,2) block being zero. Clearly, the complexity of the structure of the building is reflected in the size of the stiffness matrix, but also in the number of constraints, and thus the constraint matrix is large. Classic iterative solvers do not give satisfactory results for this problem class, which explains the need for new scalable iterative solvers for the solution of the linear systems.

Our approach to provide an efficient solver is based on the Golub-Kahan bidiagonalization (GKB) method, which has been widely used in solving least-squares problems and in the computation of the singular value decomposition of rectangular matrices [1]. In this talk, we will present how a variant of the Golub-Kahan bidiagonalization algorithm can be used as an iterative solver for indefinite matrices described above. Let now M be the positive-definite (1,1)-block of the augmented system, corresponding to the elasticity stiffness matrix. For each iteration of the generalized GKB method, a linear system Mz=b has to be solved. Using an augmented Lagrangian approach as in [2], we modify the (1,1)-block and find that, for a good choice of the parameter involved, the condition number of the generalized GKB system matrix is small and bounded above. We obtain as first desirable property that, when using a direct solver for the inner system, the algorithm converges in only a few iterations. Second, this result is particularly interesting for our industrial application, as the number of generalized GKB iterations needed to solve the augmented matrix system is independent of the problem size obtained from the finite element discretization.

For our application at hand, the matrix M is large, and thus a Cholesky decomposition to be done once at the beginning is no longer feasible. This step therefore requires an iterative solver, which gives an inner-outer iterative scheme. Its performance, in terms of computational cost, depends on the choice and accuracy of the inner iterative solver. In this talk, we will discuss the combination of stopping criteria for inner and outer iterations, such that the favourable property of requiring only a few iterations is maintained. Also, we will look at when to stop the algorithm to obtain the same order of accuracy as for having a direct inner solver, and furthermore, to achieve the same order of accuracy of the finite element discretisation.

[1] M. Arioli. Generalized Golub-Kahan bidiagonalization and stopping criteria. SIAM J. Matrix Anal. Appl., 34(2):571–592, 2013. [2] G. Golub, C. Greif. On solving block-structured indefinite linear systems. SIAM J. Sci. Comput., 24(6):2076-2092, 2003.

10:25-12:30 Session 15A: Applications. Room: Bighorn B.
10:25
An implicit approach to phase field modeling of solidification for additively manufactured materials
SPEAKER: Chris Newman

ABSTRACT. We develop a fully-coupled, fully-implicit approach to phase field modeling of solidification for additively manufactured materials. Predictive simulation of solidification in pure metals and alloys remains a significant challenge in the field of materials science, as micro-structure formation during the solidification of a material plays an important role in the properties of the solid material. Our approach consists of a finite element spatial discretization of the fully-coupled nonlinear system of partial differential equations at the microscale, which is treated implicitly in time with a preconditioned Jacobian-free Newton-Krylov (JFNK) method. The approach allows timesteps larger than those restricted by the traditional explicit CFL limit on structured and unstructured 2D and 3D meshes, is algorithmically scalable and efficient due to an effective preconditioning strategy. In addition, we discuss development of a computational tool based on this approach in which a subscale simulation can be utilized to directly inform a process simulation code at the continuum scale, or utilized to support an analytical model at the continuum scale. The computational tool leverages MPI and OpenMP, demonstrates parallel scalability and efficiency and is uniquely targeted to emerging exascale architectures.

10:50
Fixed-point iterative solution of the microlaser mode equation using a black-box Maxwell's equations solver via implicit Newton step calculation
SPEAKER: Wonseok Shin

ABSTRACT. Micrometer-sized laser sources, or microlasers, are crucial building blocks of modern integrated optical circuits. The operation of microlasers is accurately described by the steady-state ab-initio laser theory (SALT), which is a nonlinear eigenequation whose solutions are lasing modes. The SALT equation depends nonlinearly on both eigenvalues and eigenvectors. Nonlinear eigenvequations with such nonlinearity are typically solved by Newton's method.

However, the SALT equation cannot be solved directly by Newton's method, because its nonlinearity comes from the absolute values of complex eigenfunctions, which are nonholomorphic and therefore do not allow the Jacobian needed in Newton's method. The known workaround for this difficulty is to cast the eigenequation into the real domain by separating the real and imaginary parts of the equation. This approach works, but it destroys the structure of the original complex operator whose major potion is the Maxwell operator, and therefore makes it impossible to utilize the existing Maxwell's equations solvers in calculating the Newton steps.

In this work, I will present an alternative method to apply Newton's method to solving the SALT nonlinear eigenequation. The new method constructs a complex nonlinear equation (separate from the SALT equation) whose solution is equivalent to the Newton step in the real domain. This complex nonlinear equation is solved implicitly by the fixed-point iteration. It turns out that each step of the fixed-point iteration requires solving linear systems, which are in the form of Maxwell's equations and therefore can be solved by any black-box solvers Maxwell's equations solvers, including iterative solvers.

The new method has two major advantages over the previous workaround. First, compared to the previous workaround where the constructed linear system does not have any direct relationship with Maxwell's equations, the new method constructs linear systems in the form of Maxwell's equations, for which efficient iterative solution techniques (e.g., effective preconditioners) already exist. Second, compared to the previous workaround that constructs a single large linear system whose size increases with the number of laser modes participating in the laser operation, the new method constructs multiple smaller linear systems of the same size that can be solved in parallel.

The above-mentioned two advantages make the new method ideal for calculating 3D microlaser modes in parallel computing environments. The new method is implemented in the Julia programming language. Numerical examples proving the effectiveness of the new method will be presented.

11:15
Multilevel methods for the optimal control of elliptic equations with stochastic coefficients

ABSTRACT. In this work we develop multilevel preconditioners for the linear systems arising in the optimal control of elliptic equations with stochastic coefficients. The precise model problem under scrutiny is to minimize a standard regularized quadratic tracking-type cost functional averaged over the sample space, constrained by a linear elliptic equation with stochastic diffusion coefficients.

Since this model problem is a natural first candidate for an optimal control problem constrained by a stochastic partial differential equation (SPDE), it has received significant attention over the last few years. However, the majority of the research on this topic relies on solving SPDEs numerically using specialized sampling methods (multilevel Monte Carlo, sparse grid stochastic collocation); the advantage of these approaches lies in the fact that they ultimately rely on solving deterministic PDEs, albeit very many of them. An alternative approach to solving SPDEs numerically rests on using generalized polynomial chaos expansions (gPC) for modeling randomness, followed by their numerical solution via stochastic Galerkin (SG) discretizations. Such discrete formulations of SPDEs require the solution of very large, but sparse, systems of equations (linear or nonlinear). However, these systems have a structure that can be exploited to produce good solvers.

The approach to optimal control of SPDEs we adopt in this project is to formulate SPDE-constrained optimization problems based on gPC-SG discretizations. As expected, in this approach the first order optimality conditions lead to systems of equations that are even larger than those representing the discrete SPDE. In this work we show that multilevel techniques developed earlier for efficiently solving optimization problems constrained by elliptic equations can be generalized successfully to the SG formulation for optimal control of stochastic elliptic equations.

11:40
Optimal-order preconditioners for the reformulated Morse-Ingard equations
SPEAKER: Peter Coogan

ABSTRACT. The Morse-Ingard equations are a system of elliptic, complex-valued partial differential equations describing time-harmonic thermoacoustics. They are an important aspect of modeling devices to detect trace amounts of gases. We consider a finite element discretization and prove error estimates based on a Gårding-type inequality. More importantly, the finite element discretization leads to a large, sparse, ill-conditioned block system of equations. To enable scalable iterative methods, we consider block preconditioners. Both block Gauss-Seidel and Jacobi type preconditioners admit rigorous eigenvalue bounds. We also demonstrate the effectiveness with numerical results using deal.II and PETSc.

10:25-12:30 Session 15B: Numerical Linear Algebra. Room: Bighorn C/1.
10:25
Numerical methods for large-scale Lyapunov equations with symmetric banded data
SPEAKER: unknown

ABSTRACT. The numerical solution of large-scale Lyapunov matrix equations with symmetric banded data has so far received little attention in the rich literature on Lyapunov equations. We aim to contribute to this open problem by introducing two feasible solution methods, which respectively address the cases of well conditioned and ill conditioned coefficient matrices. The proposed approaches conveniently exploit the possibly hidden structure of the solution matrix so as to deliver memory and energy saving approximate solutions. Numerical experiments are reported to illustrate the potential of the described methods.

10:50
Propagation of Compression Error of ZFP with Fixed Precision for Floating-Point Data in Iterative Methods

ABSTRACT. ZFP compression is a lossy block compression algorithm that is capable of inline compression and decompression. Only the data needed for the update, at a particular value, needs to be decompressed, while the remaining data can remain in a compressed state. Since, ZFP is lossy, i.e., it allows an inexact approximation of the original data, error due to compression and decompression needs to be considered. However, numerical simulation is fundamentally about approximation, and the solution state already contains traditional errors (floating-point round-off, truncation error, and discretization error.) If ZFP’s parameters are set judiciously, lossy compression does not lose useful information if the lost bits represent other errors. If, however, a higher rate of compression is needed, the additional error caused by ZFP will contaminate the current iterate. It must be shown that these errors, as they propagate through the simulation, are bounded to prove that the process is stable. This talk first presents a static error analysis of ZFP for fixed precision. Specifically, we introduce the infinite bit vector space that allows us to work with binary representations. We can then define corresponding operators for each step of ZFP, in the infinite bit vector space, to help establish an error bound for both a static and iterative implementation of ZFP with fixed precision. To conclude, we present results from several numerical experiments to observe the error bound in a variety of settings.

11:15
Multigrid Optimization over Graph Lineages, with Applications to ANN Training
SPEAKER: Cory Scott

ABSTRACT. We introduce and develop a novel theory of process-based hierarchical model architectures, and apply it to upscaling weights between neural networks. Our work begins with the generation of graph "lineages", exponentially growing families of hierarchically related graphs. We present several methods of constructing these scale spaces, including graph generation via local graph rewrite rules and two new "skeletal" graph products we define. We characterize these model structures in relation to one another with a diffusion process- derived measure of distance which encodes similarity of two graphs with regard to a graph-local diffusion process, and explore properties of this distance. We develop methods for low-error inference of model states between members of the same model hierarchy. These methods make primary use of "prolongation matrices", maps which optimally relate graphs at coarse and fine scales. We calculate optimized prolongation maps for a variety of graph similarity problems. Finally, we demonstrate two applications of the above theory: 1) using optimized prolongations between 2D grids to generate proposal distributions in a 2D Ising (spin-glass) simulation, and 2) incorporating prolongations and their reverse operators, graph restrictions, to create a new type of "multigrid neural network" (MGNN) model that learns at multiple scales, according to classic Multi-Grid (MG) strategies of refinement and coarsification. MGNNs outperform standard autoencoder neural networks, converging to lower error in fewer training examples, by a factor of 10.