{'fraction_non_alphanumeric': 0.03578258022265881, 'fraction_numerical': 0.013503601833660773, 'mean_word_length': 4.844381506429884, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Preconditioned Least-Squares Petrov-Galerkin Reduced Order Models

This paper introduces a methodology for improving the accuracy and efficiency of reduced order models (ROMs) constructed using the least-squares Petrov-Galerkin (LSPG) projection method through the introduction of preconditioning. Unlike prior related work, which focuses on preconditioning the linear systems arising within the ROM numerical solution procedure to improve linear solver performance, our approach leverages a preconditioning matrix directly within the LSPG minimization problem. Applying preconditioning in this way can improve ROM accuracy for several reasons. First, preconditioning the LSPG formulation changes the norm defining the residual minimization, which can improve the residual-based stability constant bounding the ROM solution's error. The incorporation of a preconditioner into the LSPG formulation can have the additional effect of scaling the components of the residual being minimized, which can be beneficial for problems with disparate scales. Importantly, we demonstrate that an 'ideal preconditioned' LSPG ROM (a ROM preconditioned with the inverse of the Jacobian of its corresponding full order model, or FOM) emulates projection of the FOM solution increment onto the reduced basis, a lower bound on the ROM solution error for a given reduced basis. By designing preconditioners that approximate the Jacobian inverse, a ROM whose error approaches this lower bound can be obtained. The proposed approach is evaluated in the predictive regime on several mechanical and thermo-mechanical problems within the Albany HPC code. We demonstrate numerically that the introduction of simple Jacobi, Gauss-Seidel and ILU preconditioners into the Proper Orthogonal Decomposition/LSPG formulation reduces significantly the ROM solution error, the reduced Jacobian condition number, the number of nonlinear iterations required to reach convergence, and the wall time.


Introduction
Numerous modern-day science and engineering problems require the simulation of complex systems with tens of millions of unknowns. Despite improved algorithms and the availability of massively parallel computing resources, "high-fidelity" models are, in practice, often too computationally expensive for use in time-critical settings such as design, fast turnaround analysis and control. The situation is particularly grave in applications involving multi-query analyses (e.g., optimization, uncertainty quantification), which require simulations to be repeated many times to explore the design space or to properly characterize uncertainty. The necessary calculations can present an intractable computational burden even with the projected growth in computing power as we approach the exascale computing age.
Reduced-order modeling is a promising strategy for reducing the computational cost of such simulations while preserving high levels of fidelity. Reduced-order models (ROMs) are models constructed from high-fidelity simulations that retain the essential physics and dynamics of their corresponding full order models (FOMs), but have a much lower computational cost. Although some consider a data-fit or low-fidelity model a ROM, herein the term "reduced order model" refers to a projection-based ROM. In projection-based model reduction, the state variables are restricted to reside in a low-dimensional subspace, typically computed offline through a data-compression process performed on a set of snapshots collected from a high-fidelity simulation or physical experiment, followed by truncation. There are numerous approaches in the literature for computing a lowdimensional subspace e.g., Proper Orthogonal Decomposition (POD) [72,41], Dynamic Mode Decomposition (DMD) [62,68] balanced POD (BPOD) [61,82], balanced truncation [37,50], and the reduced basis method (RBM) [63,80]. Herein, without loss of generality, we restrict attention to the POD approach for calculating the reduced bases due to its prevalence and simplicity. Once a reduced basis is computed, the ROM dynamical system is obtained by projecting the governing equations, or some discretized form of these equations, onto the low-dimensional subspace. For nonlinear problems, an additional approximation, referred to as "hyper-reduction", is usually required to gain a computational speed-up. Since ROMs are, by construction, low-dimensional and inexpensive to evaluate, these models can enable real-time analysis and alleviate the computational burden posed by high-dimensional uncertainty quantification (UQ) problems critical to many applications.
In order to serve as a viable predictive tool, a ROM should possess certain fundamental mathematical properties, such as consistency (with respect to its corresponding high-fidelity model), stability (in space as well as in time), and convergence (to the solution of its corresponding highfidelity model). Galerkin projection, the most popular approach, is considered continuous optimal, as the resulting ROM minimizes the time-continuous residual in the l 2 norm [15], and preserves problem structure [21,22,47] for certain classes of problems (e.g., Lagrangian dynamical systems). However, the method can give rise to nonphysical instabilities [60,45,13,44,12,11,10,18] and inaccurate long-time responses [71,55,15]. Additionally, it lacks in general an a priori convergence guarantee [62]. The least-squares Petrov-Galerkin (LSPG) projection method has been proposed [17] to remedy some of these difficulties via symmetrization of the discrete Jacobian. This method performs projection at the level of the fully discrete partial differential equations (PDEs), i.e., after the PDEs have been discretized in space and time, and computes a solution that minimizes the 2 -norm of the time-discrete residual arising in each time step. This procedure ensures that adding basis vectors yields a monotonic decrease in the least-squares objective function defining the underlying minimization problem. The method can maintain efficiency through the incorporation of hyper-reduction approaches such as gappy POD [31], the discrete empirical interpolation method (DEIM) [25] or the Energy-Conserving mesh Sampling and Weighing (ECSW) method [33,24,32], introduced to maintain efficiency of the method for nonlinear problems. When combined with gappy POD, the LSPG approach is equivalent to the Gauss-Newton with approximated tensors (GNAT) method [18]. While LSPG projection does not necessarily guarantee a priori accuracy and stability, it has been shown numerically to possess better stability and accuracy properties than the Galerkin method for a variety of problems [17,18,75,77,76,15]. We note that the LSPG method has been extended in the time domain [27,57], to allow for nonlinear trial manifolds rather than linear trial subspaces [48], to operate in a domain-decomposition setting [40], and to account for physics-based constraints [67].
The aim of the present paper is to develop a methodology for improving the accuracy of ROMs constructed using the LSPG projection method for a wide range of applications through the introduction of preconditioning. As shown in [15], LSPG errors are subject to a stability constant that is dictated by the residual. The use of preconditioning within the LSPG formulation improves the conditioning of the system, which improves the stability constants. Intuitively, this aligns the 2 -norm of the residual more closely with the 2 -norm of the (time-local) state error. To motivate the main contributions of this paper, we first provide a brief overview of related work involving projection-based ROMs and preconditioning.

Overview of related work
The idea of preconditioning began to make an appearance within the model-reduction literature approximately a decade ago. The bulk of the literature on this subject falls into two categories, which we overview succinctly below. The first category of methods are aimed at developing ROM-based preconditioners for high-fidelity models. In [19,58], the authors present a novel class of ROM-based preconditioners for the iterative high-fidelity solution of transient parabolic and self-adjoint PDEs. The preconditioners are obtained by nesting appropriate projections of ROMs into the classical preconditioned Conjugate Gradient (CG) iteration within the FOM. Several authors explore the idea of developing a ROM-based preconditioner to speed-up the numerical solution of the pressure system arising in high-fidelity reservoir modeling [8,43]. In [8], Astrid et al. employ ROM concepts to develop a new preconditioner to speed up the pressure solution during a two-stage preconditioning procedure in a fully implicit time-integration scheme within a high-fidelity reservoir simulator. In [43], Jiang introduces a ROMbased preconditioner within a stationary Richardson iteration for the pressure system in a similar high-fidelity reservoir simulation.
The second category of methods explore the idea of introducing a preconditioner within the ROM workflow to speed up the online evaluation of the ROM. This work is motivated by the observation that the repeated solution of linear systems, characterized by dense matrices, can be the main computational bottleneck in efficient scaling of ROMs [69]. Commonly, the linear systems arising within the ROM workflow are solved using direct methods. This approach can be reasonable for small ROMs, but for larger ROMs containing numerous parameters, it was shown in [30] that a direct method can be more expensive to solve than the original FOM. Preconditioned iterative methods have been developed for alleviating this difficulty. In [29] and [30], preconditioners that are precomputed in the offline stage of the model reduction are introduced into a reduced basis ROM formulation with and without the application of discrete empirical interpolation to handle the nonlinear terms in the ROM. These preconditioners are derived from preconditioners used in solving the FOM, and are shown to be more efficient than direct solvers for problems depending on a moderately large number of parameters. Singh et al. develop strategies for reusing preconditioners in projection-based model reduction algorithms for parametric [70] and non-parametric [69] ROMs to improve scalability of ROM solution over direct methods. These authors introduce a novel Sparse Approximate Inverse (SPAI) preconditioner for ROMs constructed via moment matching using the Bilinear Iterative Rational Krylov Algorithm (BIRKA) algorithm. These preconditioners are aimed at systems where the linear system may change substantially between nonlinear iterations, as happens with highly nonlinear problems. Several other authors consider preconditioners constructed using IRKA methods, for which the most computationally demanding part is the solution of linear shifted systems. The ROM matrices in this case are sparse. Ahmad et al. [6] propose variants of the IRKA that use preconditioned versions of the shifted BiCG algorithm, and develop a polynomial preconditioner that can simultaneously be applied to all shifted systems. The idea of updating incomplete factorization preconditioners within the IRKA approach is explored in [7]. The use of randomized matrices to construct preconditioners for ROMs has been studied recently by several authors [83,9]. In [83], the authors propose to use randomized linear algebra for the efficient construction of a preconditioner by interpolation of matrix inverse. It is shown that, in addition to improving the online performance of a given ROM, these preconditioners can be also used to improve the quality of residual-based error estimates in the context of the reduced basis method. Building on [83], Balabanov [9] addresses the efficient construction (using random sketching) of parameter-dependent preconditioners that can be used to improve the quality of Galerkin projections or for effective error certification for reduced basis and POD ROMs with ill-conditioned operators.
It is worthwhile to remark that the idea of preconditioning an LSPG ROM is related to several other recently-introduced methodologies for improving ROM accuracy, namely the idea of scaling the ROM residual prior to its minimization and the idea of modifying the norm in which the residual is minimized. In [81], Washabaugh proposed an approach for improving the accuracy of ROMs for systems of PDEs in which different variables have drastically different scales (e.g., applications where dimensional equations are employed). Specifically, it was shown in this reference that scaling the ROM residual to get all the equations to be roughly the same order of magnitude (equivalent to the application of a diagonal preconditioner) can result in more accurate ROMs. In a similar vein, it was demonstrated in [26] that normalizing snapshot data prior to computing POD modes in order to remove dimensional effects increases numerical robustness of ROMs for hypersonic aerodynamic simulations, especially in cases with non-equilibrium thermo-chemical effects. More recently, in [56], Parish and Rizzi developed a set of dimensionally-consistent inner products and demonstrated that these inner products have a positive impact on both Galerkin and LSPG ROMs with application to the compressible Euler equations. In [28], Collins et al. examined an optimal test basis for LSPG ROMs that are obtained by solving an adjoint-like system, which transforms the residual error minimization problem to be equivalent to a least-squares state error minimization problem. The resulting test functions can be seen as incorporating the action of a preconditioner matrix that approximates the inverse of the Jacobian matrix from a ROM's corresponding FOM. In [5], Abgrall et al. proposed a residual minimization model reduction approach, in which the norm in which the residual is minimized is the L 1 norm. If the resulting L 1 residual minimization is solved iteratively using reweighted least-squares, the procedure can be interpreted as a least-squares solve.

Contributions and organization
A distinguishing feature of the present work that differentiates it from most of the references described above is the premise that the introduction of a preconditioner into the LSPG formulation can not only lead to performance gains (by reducing the condition number of the reduced Jacobian), but can actually improve a ROM's accuracy. Unlike prior work, e.g., [30,29,69], the idea herein is not to simply precondition the linear systems arising within the LSPG algorithm to improve linear solver performance, but rather to insert a preconditioner matrix directly into the minimization problem underlying the LSPG formulation. Doing this has the effect of changing the norm defining the residual minimization problem underlying the LSPG method, which can improve the residual-based stability constant bounding the method's error. We additionally demonstrate that: • an ideal preconditioned LSPG ROM (a ROM in which the preconditioner is the inverse of the Jacobian of its corresponding FOM) emulates projection of the FOM solution increment onto the reduced basis, which defines a lower bound on the error of a ROM for a given reduced basis; • it is possible to obtain a ROM whose error approaches this lower bound by improving the quality of the preconditioner employed within the LSPG formulation; • the addition of certain preconditioners into the LSPG formulation can have the effect of scaling the components of the residual being minimized to make them roughly of the same magnitude, which can minimize bias and reduce the number of nonlinear iterations required for convergence when applying the LSPG method to problems with disparate scales (e.g., dimensional PDEs, multi-physics problems); • while our approach is not specifically based on preconditioning the linear systems arising within the LSPG iteration process, it has the effect of reducing the condition numbers of an LSPG ROM's reduced Jacobians.
We evaluate the performance and predictive accuracy of several preconditioned POD/LSPG ROMs on one mechanical and two thermo-mechanical problems implemented within the open-source HPC multi-physics finite element code known as Albany [65]. We point out that the vast majority of the current ROM literature uses Galerkin projection for the reduction step of the ROM workflow when considering solids and structures [1,79,33,24,32,20,23]; herein, we demonstrate that the preconditioned LSPG method can also be effective for this class of problems. For the mechanical problem, we demonstrate that preconditioning is needed to obtain a convergent POD/LSPG ROM. For the thermo-mechanical problems, we demonstrate that, while it is possible to obtain a convergent unpreconditioned POD/LSPG ROM for certain basis dimensions, the addition of preconditioning reduces the error by between two and seven orders of magnitude, all while reducing the overall wall time by as much as 12×. The wall time improvements can be attributed to a preconditioner's ability to reduce not only the condition number of the reduced Jacobians arising within the LSPG algorithm, but also the total number of nonlinear iterations required for ROM convergence. The remainder of this paper is organized as follows. In Section 2, we provide an overview of the POD/LSPG approach to model reduction applied to a generic system of nonlinear algebraic equations, which might arise from the discretization of a PDE. We introduce the notion of preconditioning within the LSPG residual minimization problem in Section 3. We demonstrate that the use of a preconditioner has the effect of changing the norm in which the LSPG residual minimization problem is solved, and show that the use of an ideal preconditioner (the inverse of a Jacobian matrix in a particular Newton iteration) will lead to an l 2 -optimal projection of the solution increment onto the reduced basis. In Section 4, we describe succinctly the Albany multi-physics finite element code and our implementation of model reduction capabilities within this code, including a partitioned/"blocking vector" approach for applying ROM Dirichlet boundary conditions within a finite element code that does not remove the constrained (Dirichlet) degrees of freedom from the global finite element system prior to performing its numerical solution. We study the performance of the proposed preconditioned LSPG ROMs in the context of several mechanical and thermo-mechanical solid mechanics problems in Section 5, demonstrating the efficacy of the proposed preconditioning strategy. Finally, conclusions are offered in Section 6.
2 Problem formulation 2.1 Full-order model Since the methodology described herein can be applied to a wide range of FOMs arising in a variety of applications, let us consider the following generic system of nonlinear algebraic equations defining our FOM: Here, w ∈ R N is the state vector and r : R N → R N is the nonlinear residual operator. Systems of the form (1) are obtained by fully discretizing a set of governing PDEs using a spatial discretization method and, in the case of a dynamics problem, a time-integration scheme. Assuming we solve (1) using a (globalized) Newton's method, the following sequence of solutions is generated: Here, J (k) := ∂r ∂w w (k) ∈ R N ×N is the FOM Jacobian, r (k) := r w (k) ∈ R N is the FOM residual, w (0) is an initial guess for the solution, α k ∈ R is the step length (obtained using a line-search method, or set to one, as commonly done), and K is the number of Newton iterations.

Reduced order models
Our task is to build a ROM for (1) using the projection-based model reduction approach. This approach consists of three steps: (i) calculation of a reduced basis, (ii) projection of the governing equations (in our case, (1)) onto the subspace spanned by the reduced basis, and (iii) hyper-reduction to handle efficiently the projection of the nonlinear terms.
We describe each of these steps succinctly in the following subsections. As mentioned earlier in the introduction, we restrict our attention herein to reduced bases calculated using the POD. The numerical implementation of the model reduction algorithms outlined in this section within the Albany multi-physics finite element code is discussed later, in Section 5.

Reduced basis calculation via the POD
The first step in the projection-based approach to model reduction is the calculation of a basis of reduced dimension M N (where N denotes the number of degrees of freedom (dofs) in the full order model (1)) using the POD. The POD is a mathematical procedure that, given an ensemble of data and an inner product, denoted generically by (·, ·), constructs a basis for the ensemble that is optimal in the sense that it describes more energy (on average) of the ensemble in the chosen inner product than any other linear basis of the same dimension M . The ensemble {w p : p = 1, . . . , P } is typically a set of P instantaneous snapshots of a numerical solution field, collected for P values of a parameter of interest, or at P different times. Following the so-called "method of snapshots" [72], given a snapshot matrix W := [w 1 , ..., w P ] ∈ R N ×P , a POD basis Φ M of dimension M ≤ P is obtained by first computing the (thin) singular value decomposition (SVD) where the left singular vector matrix U := [u 1 , ..., u P ] ∈ R N ×P is orthogonal (U T U = I), the diagonal singular value matrix Σ = diag(σ i ) contains the ordered singular values σ 1 ≥ σ 2 ≥ ... ≥ σ P ≥ 0 of W , and the right singular vector matrix V ∈ R P ×P is also orthogonal, like U (V T V = I).
Given (3), the sought-after POD basis is obtained by selecting the first M ≤ P left singular vectors of W : It follows that Φ M has orthonormal columns and satisfies Φ T M Φ M = I. Once Φ M is calculated, we approximate the solution to (1) bỹ where φ i := u i for i = 1, ..., M ,ŵ := [ŵ 1 · · ·ŵ M ] T ∈ R M denote the generalized coordinates, and w 0 ∈ R N denotes a reference solution, often taken to be the initial condition in the case of an unsteady simulation. Substituting the approximation (5) into (1) yields (6) is a system of N equations in M unknownsŵ. As this is an over-determined system, it may not have a solution.

Least-Squares Petrov-Galerkin (LSPG) projection
As discussed earlier, the LSPG approach [17] to model reduction has shown some promise in mitigating the stability problems and other disadvantages of other approaches, e.g., the Galerkin projection method. In the LSPG approach, the ROM solution in (5) is obtained by solving the following least-squares optimization problem: Here, the approximate (Petrov-Galerkin) solution isw PG := w 0 + Φ MŵPG . The name "LSPG" ROM comes from the observation that solving (7) amounts to solving a nonlinear least-squares problem. The two most popular approaches for this are the Gauss-Newton approach and the Levenberg-Marquardt (trust-region) method. Following the work of Carlberg et al. [18], we adopt the Gauss-Newton approach 1 . This approach implies solving a sequence of linear least-squares problems of the form w (k) where K PG is the number of Gauss-Newton iterations. It can be shown that the approximation upon convergence isw PG =w where J (k) and Ψ (k) Equation (8) can be solved numerically using a variety of methods, including solving the normal equations (11), or by employing more numerically stable methods such as the QR decomposition and the SVD; see [18] for more details.
Remark 1. It is clear from (12) and (13) that LSPG projection can be interpreted as a Petrov-Galerkin process of the Newton iteration with trial basis Φ M and test basis Ψ (13).

Remark 2.
When Newton's method is employed to solve the reduced state equations for the LSPG method, the ROM solution procedure can be viewed as minimizing at each iteration the error in the search direction: where Θ := J (k) T J (k) and ||x|| Θ := √ x T Θx. Since Θ is always a symmetric positive-definite matrix, the LSPG approach has the effect of symmetrizing the discrete ROM Jacobian. As shown in [17], when (14) is satisfied, the error measure Φ M y − δw (k) Θ decreases monotonically as vectors are added to the POD basis.

Hyper-reduction
The LSPG projection approach described in Section 2.2.2 is inefficient for nonlinear problems. This is because the solution of the ROM system requires algebraic operations that scale with the dimension of the original full order model N . This problem can be circumvented through the use of hyper-reduction. A number of hyper-reduction approaches have been proposed, including DEIM [25], "best points" interpolation [53,54], collocation [49], gappy POD [31] and the ECSW method [33,24,32]. The basic idea behind these approaches is to compute the residual at some small number of points q with q N , encapsulated in a "sampling matrix" Z. This set of q points is typically referred to as the "sample mesh". The "sample mesh" is computed offline using a greedy algorithm following the notion of gappy data reconstruction [31], as described in detail in [18] and [16]. The LSPG projection approach combined with gappy POD hyper-reduction is equivalent to the GNAT method [18]. In this method, the nonlinear least-squares problem (7) is replaced witĥ where W ∈ R N ×M is a reduced basis for the residual r and the '+' symbol denotes the pseudo-inverse.
It is straightforward to demonstrate that the Gauss-Newton iterations corresponding to (15) take the form where K GNAT denotes the number of Gauss-Newton iterations. It is noted that the gappy POD approximation in (15) and (16) aims to approximate the entire residual and Jacobian (via leastsquares approximation) rather than simply sub-sample those quantities.

Preconditioned reduced order models
Having overviewed our basic ROM workflow, we now introduce the concept of preconditioned ROMs to the general LSPG formulation described in Section 2.2.2. Let M ∈ R N ×N denote a non-singular preconditioner matrix. In the present work, a preconditioned LSPG ROM is obtained by inserting M into the least-squares optimization problem (7) to yield: whereŵ PPG denotes the generalized coordinates of the preconditioned Petrov-Galerkin ROM. It is noted that M r(w) = 0 is equivalent to (1) provided M is non-singular.
Similar to the standard LSPG approach, if (19) is solved using a Gauss-Newton algorithm, one obtains the following sequence of linear least-squares problems: where K PPG is the number of Gauss-Newton iterations. In (20), the superscript (k) has been added to M (k) to indicated that, just like the Jacobian J (k) , the preconditioner may change in each Gauss-Newton iteration. We will assume that the sequence of preconditioners M (k) is non-singular. The normal equations corresponding to (20) are In the remainder of this section, we bring to light and discuss the implications of several important properties of preconditioned LSPG ROMs constructed via the approach described above. Theorem 1. An "ideal" preconditioned ROM (a ROM in which M (k) = J (k) −1 for each Gauss-Newton iteration k) emulates the projection of the FOM solution increment onto a given POD basis: Since (25) is the upper limit on the ROM accuracy given a basis Φ M , it follows that the most accurate ROM possible is realizable through the introduction of preconditioning.
Proof. Assume J (k) is full rank and let M (k) = J (k) −1 . Substituting this M (k) into (23) reduces these normal equations to Recognizing that the right-hand side of (26) contains the FOM solution increment so that, from (5), δw (k) is given by (25).
Introducing preconditioning into the LSPG formulation as in (20) has the effect of modifying the norm Θ in which the solution increment minimization is performed (14). For a generic non-singular preconditioner matrix M (k) , the LSPG ROM solution procedure can be can be viewed as minimizing at each iteration the error in the search direction (14) in the norm given by . For an ideal preconditioner (M (k) = J (k) −1 ), Θ = I, which leads to an l 2 -optimal projection of the solution increment δw (k) onto the reduced basis Φ M : Proof. It is straightforward to see that (29) gives (28).
While selecting M (k) = J (k) −1 is infeasible in practice, an important corollary of Theorem 1 is that by selecting a preconditioner matrix that approximates J (k) −1 , it is possible to improve the accuracy of a given LSPG ROM, with the ROM solution approaching the most accurate ROM solution possible (for a given basis denotes the condition number of a matrix A), it is possible to improve the condition numbers of the reduced Jacobians (24) arising in the Gauss-Newton iteration process.
We end this section by remarking that there is an additional benefit of introducing a preconditioner into the LSPG formulation (7) as in (19). It has been observed [81,26,56] that minimizing the raw (unweighted) residual r can be problematic for systems of PDEs where different variables have drastically different magnitudes, which happens frequently in various applications where dimensional equations are employed. Residual components corresponding to certain variables can be very large compared to the residual components corresponding to other variables, thereby biasing the minimization procedure. As shown in [81,26,56], scaling the ROM residual to get all the equations to be roughly the same order of magnitude can remedy this problem. From (19), it is can be seen that the introduction of a preconditioner M can have the same effect.

Implementation in the Albany multi-physics finite element code
The preconditioned LSPG ROMs described herein have been implemented with a version 3 of the Albany code base [3], described succinctly in this section. Albany is an open-source 4 C++ object-oriented, parallel, unstructured-grid, implicit finite element code for solving general PDEs, developed using the "Agile Components" code development strategy [66] with mature modular libraries from the Trilinos [39] project 5 . Over the years, Albany has hosted a number of science and engineering applications, including the Aeras global atmosphere code [73], the Albany Land-Ice (ALI) [78] ice sheet model solver, the Quantum Computer Aided Design (QCAD) [36] simulator, the ACE thermo-mechanical terrestrial model of Arctic coastal erosion [34], and the Laboratory for Computational Mechanics (LCM) [74,52,51] research code. This last project comprises Albany LCM and is specifically targeted at solid mechanics applications, such as the ones considered in this paper. A more detailed description of Albany, including a detailed description of its underlying design and the physics implemented therein, can be found in [65].
In the subsections below, we describe the main ingredients enabling the construction of preconditioned LSPG ROMs within the Albany code base. We present also a special partitioned/"blocking vector" approach for applying ROM Dirichlet boundary conditions within this code base, which does not remove the constrained dofs from the global finite element system prior to performing its numerical solution. This implementation is extensible to other finite element codes with similar boundary condition treatment.

The model reduction workflow in Albany
As described in Section 2.2, the basic workflow in building a ROM consists of three steps: calculation of a reduced basis, projection of the governing equations onto the reduced basis, and calculation of a sample mesh using hyper-reduction. The implementation of this workflow within Albany for a generic nonlinear problem is demonstrated in Algorithm 1. First, the main Albany executable is run in FOM mode to collect the relevant snapshots for a given problem, which are written to an output Exodus file. Next, the reduced basis and sample mesh are calculated in Albany as a pre-processing step by utilities known as AlbanyRBGen and AlbanyMeshReduce, respectively, using data collected from a high-fidelity simulation. The AlbanyRBGen utility employs the RBGen module within the Anasazi package of Trilinos [2] to calculate a reduced basis Φ M using the POD algorithm outlined in Section 2.2.1. While it is possible to use AlbanyRBGen to generate a scalar POD basis (i.e., a separate POD basis built for each of the unknowns to a vector-valued problem), the utility was designed to construct vector POD bases, and vector POD bases were employed in the numerical experiments in Section 5. We note that AlbanyRBGen does not orthonormalize the POD modes with respect to the mass matrix for dynamic solid mechanics problems, meaning the resulting bases are l 2 (as opposed to L 2 (Ω)) orthogonal.
Hyper-reduction in Albany is also handled as a pre-processing step using another Albany utility, AlbanyMeshSample, which creates the "sample mesh" by using either the collocation or the gappy POD approach [17,18] applied to the snapshot data (see Section 2.2.3). Following the calculation of Φ M and sample mesh, the main Albany executable is run again, this time in ROM mode, with Φ M and the sample mesh as inputs. During this step, Albany solves the optimization problem (7) (or (19) if preconditioning is used) using the Gauss-Newton method, and evaluates the ROM solution. The details of the implementation of the projection step of the model reduction procedure in Albany are described below in Section 4.2. The linear least-squares problem arising within each Gauss-Newton step is solved using the normal equations approach ((11) or (23)). Solutions to these linear systems are obtained using Krylov iterative methods (Conjugate Gradient or GMRES), implemented within the AztecOO library [38] in Trilinos.

Projection via the Albany::ModelEvaluator wrapper
The implementation of the projection step of the model reduction procedure within Albany relies heavily on the EpetraExt::ModelEvaluator nonlinear model abstraction available within this code [59], whose concrete implementation in Albany is referred to hereafter as the Albany::ModelEvaluator. The Albany::ModelEvaluator provides an extensible interface to the underlying application code (Albany) for a variety of solver and analysis algorithms within Trilinos. Suppose we are solving a nonlinear algebraic system that can be written as where r is the residual, w is the solution and µ is a vector of parameters, using the basic Newton method. The Newton solution procedure requires querying the application for various quantities, including the residual (30), and the Jacobian, defined as dr/dw, calculated in Albany using automatic differentiation within the Sacado library. The objective of the model evaluator interface is to provide these quantities to Trilinos solvers through the Albany::ModelEvaluator interface, which is agnostic to the physics, model and data structure that are needed to support the relevant matrix and vector abstractions within the solver. The Albany::ModelEvaluator has utility beyond this simple nonlinear solver example: it works in a similar fashion to perform time integration, continuation, sensitivity analysis, stability analysis, optimization and uncertainty quantification. For these more sophisticated analysis cases, it returns additional quantities besides r and dr/dw, e.g., the time-derivative of w (denoted byẇ), sensitivities of r with respect to parameters µ (denoted by dr/dµ), etc. The model reduction capabilities in Albany are implemented by creating a wrapper around the Albany::ModelEvaluator class within this code, known as the MOR::ReducedOrderModelEvaluator. This class is only activated when reduced order model options are enabled in an Albany input file, and takes the objects returned by the Albany::ModelEvaluator, e.g., the residual r (30) and Jacobian dr/dw, together with a reduced basis Φ M , and returns the reduced residual and Jacobian defined in (8) for the LSPG method. The ROM formulation within Albany is generic enough to create ROMs for various solvers and physics sets enabled in Albany; herein we focus on quasi-static mechanical and thermo-mechanical problems from the LCM suite within Albany solved using the LOCA continuation package within Trilinos [4]. There are currently two projection options implemented in Albany for the type of reduced order model, Galerkin and LSPG projection. Herein, attention is restricted to the latter approach.

Implementation of boundary conditions
In model reduction, it is customary to work with equations governing the unconstrained dofs, i.e., the dofs modulo any Dirichlet boundary conditions imposed. While in the standard finite element method formulation, constrained dofs are typically removed from the nonlinear system being solved, it is not uncommon for finite element codes to retain these dofs. This approach is taken in the boundary condition implementation within Albany. Specifically, in place of the system (1), the governing equations in Albany take the formr In (31),w ∈ RN andr : RN → RN denote the "extended" state and residual, respectively, with N = N + N DBC , where N DBC and N are the number of constrained (Dirichlet) and unconstrained dofs in the system, respectively. The extended variables can be rearranged to take the form: where w DBC ∈ R NDBC denote the constrained dofs that extend the solution w ∈ R N . Assuming (31) is solved using a globalized Newton method as in (2), Albany solves the following sequence of problemsJ (k) δw (k) = −r (k) , k = 1, . . . , K, In (34), the entries ofJ A consequence of the above discussion is that care must be taken to ensure boundary conditions are properly applied in ROMs constructed within the Albany code. When computing the reduced basis Φ M using Albany's AlbanyRBGen utility (step 2 of Algorithm 1), all constrained dofs (dofs for which Dirichlet BCs are imposed) are excluded from the snapshot matrix on which POD is performed. The resulting basis modes are then augmented with a set of modes, denoted by Φ DBC ∈ R NDBC×MDBC , defined on a so-called blocking region of the domain when performing the projection-based reduction. These modes, Φ DBC , are used as modes for the Dirichlet dofs, with common dofs grouped together into blocks. For instance, individual components (x, y, z) of a certain nodeset might be grouped into a block. Specifically, where the columns of Φ DBC ∈ R NDBC×MDBC are constructed from the (normalized) summation of the canonical unit vectors of each block of common Dirichlet dofs, i.e. Φ DBC := [u 1 , u 2 , . . . , u MDBC ], where u i := b i / b i , b i = j∈Ci e j , and C i is the set of dofs for each block of dofs comprising the blocking region as described above.
Performing an LSPG projection amounts to solving the following system where andŵ PG ∈ R M andŵ DBC ∈ R MDBC are the ROM solution (equations (5) and (7)) and blocking solution, respectively.

Preconditioned ROMs in Albany
To precondition the LSPG ROMs constructed by the Albany code base, we rely on the Ifpack suite of preconditioners [64] within Trilinos. The specific preconditioners tested herein are summarized in Table 1. In all Ifpack preconditioned cases, the effect of preconditioning is to pre-multiply the terms J (k) Φ and r (k) in (8) by an approximation of the FOM Jacobian inverse, M (k) ≈ J (k) −1 . We limit our attention to three Ifpack preconditioners: a Jacobi preconditioner, a Gauss-Seidel preconditioner and an incomplete LU (ILU) preconditioner with a level-of-fill of 1. The Jacobi preconditioner approximates J (k) −1 by an inverse of the diagonal of the Jacobian J (k) , whereas the Gauss-Seidel preconditioner is constructed by inverting the elements of the upper triangular part of J (k) . Our ILU preconditioner is another simple preconditioner, obtained by performing an incomplete LU factorization of J (k) and allowing only one non-zero diagonal above and one non-zero diagonal below the main diagonal. For a detailed description of these preconditioners and their implementations within Ifpack, the interested reader is referred to [64]. By considering several preconditioners, we are able to study numerically how the ROM error changes as the preconditioner is M (k) is improved. In addition to the preconditioners provided by Ifpack, we consider also, when computationally feasible to calculate, the ideal preconditioned case, in which M (k) is calculated by computing the action of the inverse of the Jacobian matrix J (k) −1 . As shown earlier, in Theorem 1, calculating M (k) = J (k) −1 is equivalent to directly computing the FOM solution increment δw k by solving J (k) δw k = −r (k) , then projecting that solution increment onto the ROM space to obtain the ROM solution increment δŵ (k) PG , as in (25). Since this approach solves the FOM system at every iteration instead of a modified ROM system as with the other approaches, the dominant source of error with this approach will be in the projection, i.e., how much the basis Φ M has been truncated, and we expect the errors with this approach to be very low. 5 Numerical examples

Governing equations
The predictive ability of the Albany model reduction capabilities described in Section 4 are evaluated on three problems based on the quasi-static mechanics and thermo-mechanics PDEs implemented within Albany's LCM suite. The last of these cases is a problem of practical scale and involves a realistic geometry. Before presenting our test cases, we briefly summarize the PDEs being solved.

Mechanical formulation
Consider a body defined by an open set Ω ⊂ R 3 undergoing a motion described by the mapping x = ϕ(X) : Ω × R 3 , where X ∈ Ω. Let F := Grad ϕ be the deformation gradient. Let also ρ 0 B : Ω → R 3 be the body force, with ρ 0 the mass density in the reference configuration. It is straightforward to show [52] that the Euler-Lagrange equations for a canonical quasi-static mechanical problem take the form: Div P + ρ 0 B = 0, in Ω.
Here, P = ∂W/∂F denotes the first Piola-Kirchhoff stress and, assuming an elastic material model 6 , W = W (F ) is the Helmholtz free-energy density. Assume that the boundary of the body is Γ = ∂ ϕ Ω∪∂ T Ω, where ∂ ϕ Ω is a prescribed position boundary, ∂ T Ω is a prescribed traction boundary, and ∂ ϕ Ω ∩ ∂ T Ω = ∅. Denoting the prescribed boundary positions or Dirichlet boundary conditions by χ : ∂ ϕ Ω → R 3 , and the prescribed boundary tractions or Neumann boundary conditions are T : ∂ T Ω → R 3 , (38) is supplemented by the following constraints: Embedded within the Helmholtz free-energy density W (F ) is a constitutive or material model. For the numerical examples given here, we employ a Neohookean constitutive model extended to the compressible regime, which is consistent with a nonlinear elastic material [42]. Following common practice in solid mechanics, we decompose W (F ) into its volumetric and deviatoric components: For a Neohookean material, the components comprising W (F ) in (40) are given by: In (41), κ denotes the bulk modulus, given by where E is the elastic modulus and ν is Poisson's ration; µ denotes the shear modulus, given by J := det(F ) is the determinant of the deformation gradient, and is the left Cauchy-Green deformation, tensor, where F M denotes the mechanical part of the deformation gradient. For a pure mechanics problem, F M = F . The mechanical equations (38) are solved for the displacements, the primary unknowns in these equations.

Thermo-mechanical formulation
In the case of a quasi-static thermo-mechanical problem, the mechanics equations (38) are augmented with the following generic steady-state heat conduction equation: where q denotes the heat flux, T denotes the temperature and f is a source term. Equation (45) can include boundary conditions of both the Dirichlet and Neumann type, e.g., where T p and q p are prescribed temperature and heat flux values, respectively. In this case, a temperature dependence is introduced into the Helmholtz free-energy density (41): W = W (F , T ). For the Neohookean formulation considered herein, temperature effects are introduced in the form of thermal expansion. More specifically, the mechanical part of the deformation gradient appearing in (44) now takes the form: where is the thermal part of the deformation gradient, with I denoting the identity tensor. In (48), α is the thermal expansion coefficient and T ref is the reference temperature. For the problems considered here, the heat flux in (45) takes the form: where K is the thermal diffusivity tensor.

Evaluation criteria
Recall that we denote the FOM solution as w and the approximate ROM solution asw = w 0 +Φ Mŵ (5). Let w i andw i denote the FOM and ROM solutions, respectively, corresponding to time t i , for i = 1, ..., P , where P denotes the number of snapshots collected. Building on these definitions, we define the global relative error in the approximate ROM solution as: For mechanical problems, w andw represent the displacement field solutions, whereas for thermomechanical problems, these vectors include also the temperature field. In the subsequent sections, the proposed preconditioned LSPG ROMs are evaluated using the metric (50). We emphasize that the evaluations performed herein are in the predictive regime. That is, we calculate errors by running the proposed ROMs with different parameter values than the parameter values used to generate the training data from which the ROMs were constructed. More details on the training/testing for each problem considered are provided in the subsections below. Following the usual convention in solid mechanics, the FOMs and ROMs for all three test cases considered are run dimensionally. In addition to reporting the global relative error for each of the ROMs evaluated, we also report the total CPU-time of the online stage of the model reduction procedure for each ROM. The specifications of the run-time environment for each benchmark problem are described in the following subsections. As the purpose of this work is to present a preliminary study examining the viability of preconditioning within the LSPG framework, hyper-reduction, which would potentially introduce another source of error, was not employed in our study. Evaluating the performance of LSPG ROMs with preconditioning and hyper-reduction will be the subject of future work. It is noted that special hyper-reduction techniques, e.g., the ECSW method [33,24,32] of Chapman, Farhat et al. and the structure-preserving/"matrix gappy POD" method of Carlberg et al. [20,23], are needed to preserve the Lagrangian structure of ROMs for solid mechanics problems such as the ones discussed herein. These are not available within the Albany code base at the present time.
Lastly, for each problem considered, we report also either the average condition number of all reduced Jacobians obtained during the Gauss-Newton iteration process or the total number of nonlinear iterations required for ROM convergence. The latter metric is adopted for the larger thermo-mechanical pressure vessel problem (Section 5.4), as estimating reduced Jacobian condition numbers is not feasible for a problem of that size.

Mechanical and thermo-mechanical beam
Our first example involves a simple three-dimensional (3D) beam geometry Ω of size 0.16 × 0.016 × 0.032 m centered at (0, 0, 0), discretized using 320 hexahedral elements, which gives rise to 525 nodes ( Figure 1). The geometry Ω is divided into five material blocks, as shown in in Figure 1(a). The purpose of having several material blocks is to enable the specification of different material models and/or material parameters in different parts of the domain. As mentioned earlier, we consider a nonlinear elastic material model known as the Neohookean model, which is prescribed in all the blocks comprising Ω; however, as discussed below, different material parameters are specified in different blocks.
In order to complete the problem definition, it is necessary to specify well-posed boundary conditions on our beam geometry Ω. We prescribe boundary conditions on four of Ω's boundaries, shown in Figure 1(b) and denoted by Γ w , Γ b , Γ t and Γ n . Γ w corresponds to x = −0.08 and is shown in red; Γ b corresponds to y = −0.008 and is shown in yellow; Γ n corresponds to z = −0.016 and is shown in blue; Γ t corresponds to y = 0.008 and is shown in green. In both versions of our problem formulation (mechanical and thermo-mechanical; Sections 5.3.1 and 5.3.2, respectively), Γ w is fixed in the x-direction, Γ n is fixed in the z-direction and Γ t is fixed in the y-direction. Additionally, the following linearly-varying time-dependent pressure Neumann boundary condition is applied on Γ b P (t) = (7.2599 × 10 4 )t Pa, where t is a pseudo-time variable, described in more detail below. For the thermo-mechanical version (a) Discretization of Ω (b) Ω's boundaries of the beam problem (Section 5.3.2), an additional temperature boundary condition is prescribed on Γ w ; we defer discussion of this boundary condition until Section 5.3.2. For both variants of this problem, the mechanical and the thermo-mechanical beam, the system, initially at rest, is solved quasi-statically by performing a homotopy continuation with respect to the pseudo-time variable t appearing in (51). This amounts to incrementing the applied boundary conditions (e.g., the pressure load (51)) with respect to t.

Mechanical beam
We will first investigate the purely mechanical response of the beam described above; hence, the governing PDEs are those described in Section 5.1.1. For this problem, the displacement boundary conditions represent 235 constrained dofs out of the 1575 total, leaving this problem with 1340 free dofs.
For the study performed herein, different material parameters are specified within different sets of blocks. Let block B a denote the union of blocks 1, 2 and 4 (shown in green, yellow and cyan, respectively, in Figure 1(a)) and let block B b denote the union of blocks 3 and 5 (shown in magenta and blue, respectively, in Figure 1(a)). In block B a , we set the Young's modulus, Poisson's ratio and density to the following values: E a = 1.103 × 10 11 Pa, ν a = 0.32 and ρ a = 7920 kg/m 3 , respectively. The values of these same parameters in block B b , denoted by E b , ν b and ρ b , were varied for the purpose of ROM training and testing, are provided in Table 2. These values were generated by performing Latin Hypercube (LHC) sampling within the following parameter ranges:  Table 2. It is emphasized that the parameter variations in Table 2 led to nontrivial differences in the displacement field, which varied by as much as 20% between the different training and testing cases. The main results for the mechanical beam problem are summarized in Figures 2-5. The reader can observe that the classical (unpreconditioned) LSPG solution does not appear in these plots. The LSPG solution is not included in our results summary because the unpreconditioned LSPG ROMs were not able to converge to a nonlinear solver tolerance smaller than O(1) for any of the basis sizes considered. Figure 2 plots the ROM global relative error (calculated using (50)) as a function of the basis dimension for various preconditioned LSPG ROMs, including the projected solution increment ROM, indicated with a black dashed line. It can be seen that, in general, as the preconditioner is improved (from Jacobi, to Gauss-Seidel, to ILU), the accuracy of the resulting ROM solution is improved, approaching the ideal projected solution increment solution, as expected from the discussion in Section 3. Next, in Figure 3, we report wall times for the four ROMs evaluated for each of our testing cases. The wall times reported are averaged over four processors of a Linux workstation having 20 Intel Xeon CPU E5-2670 v2 CPUs. The reader can observe that the more sophisticated ILU preconditioner results in slightly larger wall times in general compared to the Jacobi and Gauss-Seidel preconditioners. This is due to a larger preconditioner construction time associated with the ILU preconditioner. As expected, the ideal preconditioner solution, which requires calculating the action of J (k) −1 in every iteration as discussed in Section 4.4, is in general the costliest to obtain. The results in Figures 2 and 3 are combined in Figure 4, which shows a Pareto plot for the mechanical beam problem, in which the global relative error is plotted as a function of the wall time. From this plot, it is possible to identify the optimal preconditioner to use, based on one's error and CPU-time requirements. The blue line in this plot traces points having Pareto optimality, which define the so-called Pareto front. Finally, in Figure 5, we report the average condition number of the reduced Jacobian, J (k) PPG (see (23)) over all Gauss-Newton iterations for each of the testing cases summarized in Table 2. Although the LSPG ROMs were not convergent for this problem, we include condition numbers for the LSPG ROMs in this plot. The reader can observe that, as expected, the addition of a preconditioner M (k) improves the condition number of the corresponding ROM relative to the baseline LSPG ROM. For all ROMs evaluated, condition numbers grow steadily with the basis dimension M ; however, preconditioning using the relatively simple preconditioners considered herein (the first three rows of Table 1) is able to reduce the condition number of the reduced Jacobian of the LSPG ROM by up to an order of magnitude. By design, the ideal preconditioner (equivalent to the projected solution increment) gives rise to a ROM system having a perfect condition number of one for all basis dimensions.  Table 2.

Thermo-mechanical beam
We now turn our attention to the thermo-mechanical version of the beam problem, in which the mechanical equations of Section 5.1.1 are augmented with the thermal equation in Section 5.1.2. In the current Albany implementation, the thermal and mechanical equations are coupled monolithically within the code, with the coupling occurring at the level of the material model, as described earlier in Section 5.1.2.
In addition to the displacement and pressure (51) boundary conditions described previously in Section 5.3.1, we impose a temperature Dirichlet boundary condition on Γ w in Figure 1(a). The value of the prescribed temperature varies linearly from 293 to 393 K over the course of the simulation: As before, the variable t in (52) is a pseudo-time variable which is incremented quasi-statically via  Table 2.
homotopy continuation to emulate time-dependent behavior in a quasi-static framework. Note that the temperature boundary condition (52) is not the only one advanced in time quasi-statically; the pressure Neumann boundary condition (51) is also incremented, as discussed earlier in Section 5.3.1. As before, we use the Neohookean material model for the mechanical problem, and specify different material properties in block B a and B b . In block B b , the mechanical properties are E a = 1.103 × 10 9 Pa (Young's modulus), ν a = 0.32 (Poisson's ratio) and ρ a = 7.92 × 10 −5 kg/m 3 (density). As discussed in Section 5.1.2, the Helmholtz free-energy density A appearing in (38) is a function of the temperature T in the thermo-mechanical case. This thermo-mechanical formulation requires one to specify a reference temperature. In block B a , we use a reference temperature of T a,ref = 293 K. In block B b , the mechanical properties of the Neohookean material as well as the reference temperature, denoted by T b,ref are varied for both the training and the testing cases, as summarized in Table 3. Additionally, we assume the material is isotropic, i.e., K := kI 3 in (49),  Table 2 are used to generate this figure.
where I 3 is the 3 × 3 identity matrix, with thermal diffusivity k a = k b = 1.1870 × 10 3 m 2 /s in both blocks. We also assume the thermal expansion coefficient (48) is the same in both blocks, taking the value α a = α b = 1 × 10 −5 K −1 . There are no applied forces and the system is initially at rest, with an initial temperature of 293 K. For the quasi-static advancement of the system, we use a step size of 1 s and run the problem for a total of 720 steps. The displacement and temperature boundary conditions represent 260 constrained dofs out of the 2100 total dofs, leaving this problem with 1840 free dofs. As for the mechanical version of this problem, the parameter ranges in Table 3 are selected such that nontrivial variations in the solution are observed; specifically, the displacement varies by as much as 60% as the material properties are perturbed according to Table 3.
The ROM was trained with six sets of parameters in block B b , chosen using LHC sampling and summarized in Table 3. The following parameter ranges for the LHC sampling were employed: . This training generated a total of 4318 snapshots, which were used to build various LSPG/POD ROMs (with and without preconditioning) ranging in size from 1 to 721 modes, like in the mechanical variant of this problem. These ROMs were tested in the predictive regime over four additional sets of parameters, also summarized in Table 3.  Table 2. Table 3: Thermo-mechanical beam problem: summary of parameters specified in the material model in block B b for the training and testing stages of the ROM process. The main results for the thermo-mechanical beam problem are summarized in Figures 6-10. It is interesting to observe that, despite the multi-physics and multi-scale nature of this problem, the classical (unpreconditioned) LSPG ROMs, plotted in red in these figures, deliver convergent solutions for the smaller basis dimensions considered. Figures 6 and 7 report the global relative error in each of the ROM solutions, computed using (50). The reader can observe by examining this figure that, while convergence with basis refinement is observed up to a point (until M ≈ 7 modes) for the unpreconditioned LSPG ROM, for larger basis dimensions, the error grows steadily before reaching a point where a lack of convergence is observed. The results are markedly different for the preconditioned ROMs. By introducing preconditioning, it is possible to reduce by between two and six orders of magnitude, depending on the basis dimension, and there are no convergence issues, like for the unpreconditioned cases. Moreover, one can see that all the preconditioned LSPG ROMs achieve errors which are close to (less than one order of magnitude greater than) the error obtained by the projected solution increment ROM, which represents an ideal preconditioned LSPG ROM.
We next investigate the effect of preconditioning on the overall run time of the ROM problem ( Figure 8). As for the mechanical version of this problem, we report wall times averaged over eight processors of a Sandia National Laboratories' Linux cluster known as Uno, which has 201 dual-socket eight-core 2.7 GHx Intel Sandy Bridge CPUs. The reader can observe from Figure 8 that for the majority of testing cases and basis dimensions, the preconditioned LSPG ROMs considered achieve smaller wall time than the baseline (unpreconditioned) LSPG ROM. Combining these results with our earlier error results (Figures 6-7), we can confidently conclude that preconditioning is extremely advantageous for the thermo-mechanical beam problem: in general, for the same computational cost, one is able to achieve a much smaller global relative error with preconditioning. This result is confirmed in the Pareto plot shown in Figure 9. As before, the black line in Figure 8 represents the ideal preconditioned ROM. This ROM has a much higher computational cost than the other ROMs evaluated, as expected.
Finally, in Figure 10, we report average condition numbers of the reduced Jacobians J (k) PPG (see eq. (23)) encountered during the Gauss-Newton iteration process for each of the ROMs considered. The most striking observation that can be made from this figure is that the regular LSPG Jacobians are very ill-conditioned, with condition numbers ranging between O(10 14 ) and O(10 20 ) depending on M , the basis dimension. These condition numbers are between seven to ten times greater than the condition numbers for the mechanical variant of this problem (see Figure 5). The extreme ill-conditioning exhibited by the thermo-mechanical beam problem can be attributed to the extreme differences in scales between the thermal and the mechanical problems: whereas the temperature solution is O(100), the displacement solution varies between O(10 −7 ) and O(10 −5 ) during the duration of the simulation, a difference of between seven and nine orders of magnitude. Our results demonstrate that, by introducing a simple preconditioning strategy into the LSPG formulation, it is possible to bring down the reduced Jacobian condition numbers by as many as ten orders of magnitude, thereby alleviating to a large extent the scaling issue introduced by the large difference in scales between the thermal and mechanical problems (a similar result was discovered by Washabaugh in [81], but for a different application, namely computational fluid dynamics). Experience involving other solid mechanics applications within Albany [52] suggests that further improvements can be made by introducing physics-specific block preconditioners into the framework; however, a study of this sort goes beyond the scope of this work. As before, the condition number of the reduced Jacobian for the projected solution increment ROM is identically one regardless of the basis dimension.
(a) Testing case 1 (b) Testing case 2 Figure 6: Thermo-mechanical beam problem: global relative errors (50) for various ROMs as a function of the basis dimension for testing cases 1-2 in Table 3.

Thermo-mechanical pressure vessel
The last test case considered aims at studying the thermo-mechanical response of a 3D pressured cylindrical vessel shown in Figure 12(b). This problem is substantially larger and more realistic than the beam problem considered earlier in Section 5.3. Our approach herein is to model a quarter of the pressure vessel (denoted by Ω), as shown in Figure 12(a), and apply appropriate symmetry boundary conditions on the relevant sides to emulate a simulation on the pressure vessel in its entirety. We   Table 3. discretize our domain Ω using a mesh comprised of 77,768 hexahedral elements, which give rise to 92,767 nodes. Similar to the beam problem, the domain Ω is comprised of five element blocks, depicted in green (block 1), yellow (block 2), blue (block 3), magenta (block 4) and cyan (block 5), as shown in Figure 11, so as to enable to specification of different materials and/or material properties in different parts of the domain. Figure 12 depicts the boundaries on which boundary conditions are applied for the thermo-  Table 2.
mechanical pressure vessel problem. Let Γ denote the boundary of our computational domain Ω. As shown in Figure 12(a), we decompose Γ as Γ := Γ 1 ∪ Γ 2 ∪ Γ 3 ∪ Γ 4 ∪ Γ 5 where Γ i ∩ Γ j = ∅ for i, j = 1, ..., 5 and i = j. In Figure 12(a), Γ 1 is shown in green, Γ 2 is shown in blue, Γ 3 is shown in yellow, Γ 4 is shown in gray and Γ 5 is shown in red. It is noted that Γ 5 consists of a single point located in the center of the top of the pressure vessel. Since we are modeling a quarter of the pressure vessel, symmetry boundary conditions are applied on the displacements on Γ 1 , Γ 2 and Γ 5 . These BCs emulate performing a simulation on the full cylindrical pressure vessel geometry ( Figure  12(b)). In addition to the displacement boundary conditions mentioned previously, we impose the time-dependent pressure Neumann boundary condition (51) and a time-dependent temperature Dirichlet boundary condition (52) on Γ 3 . This amounts to heating and pressurizing our vessel from the inside. The boundary conditions account for a total of 5020 constrained dofs out of the 371,068 total dofs, leaving this problem with 366,048 free dofs.  Table 3 are used to generate this figure.
Similar to the beam problem (Section 5.3), we group together different element blocks for the purpose of specifying the material model parameters for the mechanical problem's material model, which we take to be Neohookean. Let B a denote the union of blocks 3 and 5 (magenta and blue, respectively, in Figure 12) and let B b denote the the union of blocks 1, 2 and 4 (green, yellow and cyan, respectively, in Figure 12). We fix the mechanical parameters in block B a to have the following values for the Young's modulus, Poisson's ratio, and density, respectively: E a = 1.103 × 10 9 Pa, ν a = 0.32, and density ρ a = 7.92 × 10 −3 kg/m 3 . The thermal equation requires us to specify a reference temperature, which is set to T a,ref = 293 K within the thermal model. Similar to the thermo-mechanical version of the beam problem (Section 5.3.2), these analogous parameters in block B b are varied using LHC for the purpose of training and testing our ROMs; the values of these parameters are given in in Table 4 and were sampled from the following ranges: The following additional parameters are also specified within the thermal model in both blocks B a and B b . As for the thermo-mechanical beam problem (Section 5.3.2), we assume the materials are isotropic in both blocks (K := kI 3 ), with thermal diffusivity k a = k b = 1.1870 × 10 3 m 2 /s. Additionally, we assume the thermal expansion coefficient (48) is the same in both blocks, and prescribe it the value α a = α b = 1 × 10 −5 K −1 . There are no applied forces, and the system is initially at rest with an initial temperature of 293 K. The system is advanced using a quasi-static approach, in which the a homotopy continuation is performed with respect to the pseudo-time variable t in (52). The simulation proceeds for a total of 720 steps with a step size of 1 s.
In this example, we trained the ROM over four sets of parameters, summarized in the top part of Table 4. Snapshots were collected every continuation step, and POD bases of various numbers of modes, from 2 to 790, were constructed using our ensemble of 2160 snapshots. Once our ROMs  Table 3.
were constructed, they were evaluated in the predictive regime on two additional cases, summarized in the bottom part of Table 4. The final displacement FOM solutions and configurations for the two testing cases are plotted in Figure 13. The reader can observe that the solutions are noticeably different, indicating that our parameter variations led to nontrivial changes in the problem solution.      Figure 13: Thermo-mechanical pressure vessel problem: final displacement and configuration for the two testing cases in Table 4, obtained by running Albany in FOM mode.
The main results for thermo-mechanical pressure vessel problem are summarized in Figures 14-18. Similar to the thermo-mechanical beam problem considered in Section 5.3.2, the unpreconditioned LSPG ROMs fail to converge for the two larger basis dimensions evaluated. The ideal preconditioned ROM solution is not plotted, as calculating this solution is too computationally expensive for this problem. Figure 14 plots the global relative errors (50) as a function of the basis dimension for the two testing cases of interest. The reader can observe that the unpreconditioned LSPG approach is unable to deliver a solution with a relative error smaller than O(10%). In contrast, all three preconditioned LSPG ROMs produce solutions that are between four and six orders of magnitude more accurate. This result is confirmed by Figure 15, which plots the final displacement and configuration obtained using the FOM and two LSPG ROMs having 32 modes, an unpreconditioned LSPG ROM and an LSPG ROM preconditioned using the Gauss-Seidel preconditioner, for testing case 1. It can be seen that the preconditioned LSPG ROM (Figure 15(c)) is indistinguishable from the FOM (Figure 15(a)); in contrast, the unpreconditioned ROM solution (Figure 15(b)) is visibly incorrect. Turning our attention back to Figure 14, it can be seen that, by improving the preconditioner, it is possible to improve solution accuracy. Figure 16 depicts the wall times required to run the various LSPG ROMs evaluated for each of the testing cases. Our simulations for the thermo-mechanical pressure vessel problem were performed on 64 cores (4 nodes) of the Skybridge high-capacity cluster located at Sandia National Laboratories, which contains 1848 nodes, each having 16 2.6 GHz Intel Sandy Bridge processors. It can be seen from this figure that a speedup of as large as 12× is attainable through the introduction of preconditioning into the LSPG ROM formulation. Figures 14 and 16 together demonstrate that LSPG ROMs are substantially faster and more accurate than unpreconditioned LSPG ROMs. This result is reinforced by Figure 17, which depicts a Pareto plot and front for the thermo-mechanical pressure vessel problem. As before, it is possible to use this figure to determine the ideal preconditioner to use based on one's error and CPU-time requirements.
Our final figure, Figure 18, plots the total number of nonlinear iterations required to attain a convergent solution for the various ROMs evaluated in this section. We present this result in lieu of the average condition number of the reduced Jacobian J (k) PPG , as the thermo-mechanical pressure vessel problem is too large to calculate the latter quantity. Comparing Figure 18 with Figure 16,  Table 4. (a) Testing case 1 (b) Testing case 2 Figure 16: Thermo-mechanical pressure vessel problem: wall times (in s, averaged over 64 processors) for various ROMs as a function of the basis dimension for the two testing cases in Table 4.  Table 3 are used to generate this figure.
one can see that the wall time improvements achieved through the introduction of preconditioning are likely due to the reduction in the total number of nonlinear iterations required to achieve a convergent ROM solution: by preconditioning the LSPG formulation, one can reduce the number of nonlinear iterations by up to a factor of 12×. It is interesting to observe that, unlike for the classical LSPG formulation, the number of nonlinear iterations remains approximately constant regardless of the basis dimension for the preconditioned LSPG ROMs. The reduction in the number of nonlinear iterations when preconditioning is employed can be attributed to the scaling effect of the ROM preconditioner for this multi-scale and multi-physics problem.
(a) Testing case 1 (b) Testing case 2 Figure 18: Thermo-mechanical pressure vessel problem: total number of nonlinear iterations for various ROMs as a function of the basis dimension for the two testing cases in Table 4.

Conclusions
Projection-based model reduction can enable real-time and multi-query analyses in a variety of applications. The LSPG method to building projection-based ROMs has shown particular promise, as it can generate stable and accurate solutions for applications where standard Galerkin techniques have failed [17,18,75,77,76,15]. Despite its superior performance over Galerkin projection, we have observed that, when run in the predictive regime, the LSPG method sometimes exhibits a lack of convergence or delivers a solution whose error surpasses the engineering tolerance required by the application. The approach is most likely to suffer from these deficiencies when applied to problems with disparate scales, such as dimensional PDEs or multi-physics problems [81,26,56]. This is because, when there are differences in scale between components of a PDE system, residual components corresponding to certain variables can be very large compared to the residual components corresponding to other variables, which can bias the residual minimization procedure underlying the LSPG formulation.
In this paper, we have demonstrated that the accuracy, robustness, efficiency and convergence properties of the LSPG method for projection-based model reduction can be improved substantially through the introduction of preconditioning. While this approach can reduce the condition number of the matrix problems arising in the LSPG iteration procedure, it is unlike other existing ROM preconditioning approaches (e.g., [30,29,69]), as the preconditioning is introduced directly into the LSPG residual minimization problem. Doing so ensures that all residual components are on approximately the same scale, even when the method is applied to multi-scale and multi-physics problems having variables of drastically different magnitudes. We demonstrated that preconditioning the LSPG method in this way is equivalent to modifying the norm in which residual minimization is performed, which can improve the residual-based stability constant bounding the method's error. We additionally showed that, by designing a preconditioner that approximates the inverse of the Jacobian matrix, it is possible to create an LSPG ROM whose accuracy approaches the upper limit on ROM accuracy for a given reduced basis: an ideal preconditioned ROM, which is equivalent to the projection of the FOM solution increment onto the reduced basis.
We studied the efficacy of the proposed approach using one mechanical and two thermo-mechanical examples simulated using the Albany HPC multi-physics finite element code [65], into which we have implemented the POD/LSPG method to model reduction. Our test cases include a realistic simulation involving a thermo-mechanical vessel pressurized from the inside, having approximately 370,000 dofs. We developed a partitioned/"blocking vector" approach for applying ROM Dirichlet boundary conditions strongly within this code, which does not remove the constrained dofs from the global finite element system prior to performing its numerical solution. We considered three simple preconditioners, a Jacobi preconditioner, a Gauss-Seidel preconditioner, and an ILU preconditioner, implemented within the Ifpack Trilinos library, on which Albany has a dependence. Our ROMs were evaluated in the predictive regime, with prediction performed across the material parameter space. For the first test case, known as the mechanical beam problem, the classical LSPG method failed to deliver a convergent solution regardless of the basis dimension. In contrast, all three preconditioned ROMs converged robustly and efficiently, with global relative errors as small as O(0.01%). While the classical LSPG method achieved convergence for the smaller basis dimensions considered in our two thermo-mechanical test cases, we demonstrated that the introduction of preconditioning can reduce the global relative error by up to six orders of magnitude and the wall time by as much as 12×. The wall time gains are attributed to a reduction in the reduced Jacobian condition number and a reduction in the total number of nonlinear iterations required for convergence, both achieved via the proposed preconditioning strategy. For thermo-mechanical problems, which exhibit extreme differences in scale between the displacement and temperature solutions, the introduction of even a simple preconditioner such as Jacobi reduced the reduced Jacobian condition number by as many as ten orders of magnitude.
The results discussed herein pave the way for several subsequent studies, which will be the subject of future work. First, we plan to extend our preconditioned LSPG formulation to include hyperreduction using Lagrangian structure-preserving methods such as [33,24,32,20,23], and intend to repeat the present numerical experiments with hyper-reduction enabled. It would additionally be interesting to study: (1) the efficacy of more sophisticated preconditioners within our framework, e.g., multi-grid preconditioners that take advantage of rigid body mode information and block preconditioners for multi-physics applications, and (2) the performance of our framework with more sophisticated reduced bases, e.g., nonlinear manifold reduced bases constructed using an auto-encoder or some other type of neural network [35,48,46]. Lastly, we plan to investigate the viability of our approach for a larger set of applications, including problems in the field of Computational Fluid Dynamics (CFD) and problems in which the ROM is used to perform prediction in time. Our preliminary results [76,14] suggest that preconditioning improves accuracy and efficiency in problems involving compressible cavity flows, hypersonic aerodynamics, thermal ablation and reacting hypersonic flows.