Reduction of Nonlinear Embedded Boundary Models for Problems with Evolving Interfaces
Abstract
Embedded boundary methods alleviate many computational challenges, including those associated with meshing complex geometries and solving problems with evolving domains and interfaces. Developing model reduction methods for computational frameworks based on such methods seems however to be challenging. Indeed, most popular model reduction techniques are projectionbased, and rely on basis functions obtained from the compression of simulation snapshots. In a traditional interfacefitted computational framework, the computation of such basis functions is straightforward, primarily because the computational domain does not contain in this case a fictitious region. This is not the case however for an embedded computational framework because the computational domain typically contains in this case both real and ghost regions whose definitions complicate the collection and compression of simulation snapshots. The problem is exacerbated when the interface separating both regions evolves in time. This paper addresses this issue by formulating the snapshot compression problem as a weighted lowrank approximation problem where the binary weighting identifies the evolving component of the individual simulation snapshots. The proposed approach is application independent and therefore comprehensive. It is successfully demonstrated for the model reduction of several twodimensional, vortexdominated, fluidstructure interaction problems.
keywords:
data compression, embedded boundary method, evolving domains, immersed boundary method, interfaces, model reduction, singular value decomposition1 Introduction
Two common computational frameworks for the solution of problems with evolving domains and interfaces are the Arbitrary Lagrangian Eulerian (ALE) Donea (); CF1 (); CF2 (); ALERezoning () and Embedded Boundary Method (EBM) frameworks. In an ALE framework, the computations are performed on an interfacefitted mesh which is deformed using a mesh motion algorithm MM1 (); MM2 (); MM3 () to follow the evolution of the interface. In an EBM framework, the interface is embedded in a background mesh  also called an embedding mesh  and allowed to intersect it as it evolves (Peskin_AN_2002, ; Mittal_ARFM_2005, ; Wang_IJNMF_2011, ; FIVER, ). In the latter case, the interface separates the computational domain in two regions: a ‘‘ghost’’ (or fictitious) region usually associated with the volume delimited by the embedded interface (boundary, or body), and a ‘‘real’’ region usually associated with the remainder of the computational domain. Both computational frameworks have advantages and shortcomings. An ALE framework is particularly advantageous for problems where the evolution of the interface is characterized by smallamplitude motions and deformations. An EBM framework is often indispensable for problems where the interface undergoes largeamplitude motions and deformations, and/or topological changes.
Irrespectively of the chosen framework however, the computational cost of highfidelity simulations of problems with evolving domains and interfaces can be prohibitive. For this reason, interest in Model Order Reduction (MOR) methods continues to grow (Rewienski_LAA_2006, ; Thuan1, ; Amsallem1, ; Astrid_AC_2008, ; Chaturantabut_JSC_2010, ; Carlberg1, ; Balajewicz_JFM_2013, ; Carlberg_JCP_2013, ; Navon1, ; Navon2, ). Most of these methods are projectionbased. Therefore, they rely on the computation of suitable basis functions. Typically, these are constructed by computing solution snapshots of some instances of the family of problems of interest and compressing them using the Proper Orthogonal Decomposition (POD) or Singular Value Decomposition (SVD) method. In an interfacefitted framework, this computation is straightforward because the computational domain does not contain a ghost region. In an EBM framework however, this computation is complicated by the presence of a ghost region where the solution snapshots may or may not be populated, and in the latter case, the populated values may or may not have a physical meaning. The situation is further exacerbated when the interface evolves in time, as in this case the partitioning into meaningful and meaningless components of the collected snapshots also evolves in time, which complicates further the compression of the collected snapshots.
Therefore, the main objective of this paper is to propose an effective method for the compression of solution snaphsots collected during the simulation of problems with evolving domains and interfaces using an EBM, in view of enabling the projectionbased reduction of highfidelity embedded boundary computational models. This method consists of formulating the data compression problem as a weighted lowrank matrix approximation problem, where the weighting is binary and identifies the ghost/real partition of the individual snapshots. Therefore, the basis functions derived using this approach are optimal for the real component of the partition of the computational domain.
The weighted lowrank approximation problem formulated in this paper is pervasive. It can be found, for example, in:

Computer vision applications with occlusion.

Signal processing applications with irregular measurements in time and space.

Control problems with malfunctioning measuring devices.
This problem is a generalization of the popular compressed sensing problem (Donoho_2006, ; Candes_2008, ). It also shares many similarities with the matrix completion problem (Candes_2009, ; Cai_2010, ; Jain_ACM_2013, ; Vandereycken_SIAM_2013, ). For all these reasons, many algorithms for the solution of this problem are readily available (for example, see (Srebro_ICML_2003, ; Buchanan_2005, ; Chen_2008, )).
Whereas the approach outlined above for compressing solution snaphots computed by an EBM in view of constructing a reducedorder basis suitable for model reduction is application independent and therefore quite general, it is developed and demonstrated in the remainder for this paper for FluidStructure Interaction (FSI) problems. To this effect, the remainder of this paper is organized as follows.
In Section 2, the main notation used in this paper is laid out to facilitate its reading. In Section 3, a general EBM framework for Computational Fluid Dynamics (CFD) is overviewed, and the standard projectionbased MOR approach is recapitulated. In Section 4, the proposed method for compressing CFD snapshots computed by an EBM is introduced. In Section 5, this method is applied to three different FSI problems. Finally in Section 6, conclusions are offered and prospects for future work are summarized.
2 Notations and definitions
Throughout this paper, matrices are denoted by bold capitals (ex. ), vectors by bold lower case (ex. ), and subscripts identify rows and columns (for example, is the entry at the th row and th column of matrix ). The ‘‘colon’’ notation is used to specify columns or rows of a matrix (for example denotes the th column and denotes the first rows of matrix ). and denote the identity and null matrices, respectively.
For two matrices and of equal dimension , the Hadamard product is the matrix of the same dimension whose elements are given by
(1) 
If is an matrix and is a matrix, the Kronecker matrix product is the block matrix
(2) 
The operator denotes the vectorization of the matrix formed by stacking the columns of into a single column vector
(3) 
For a vector of dimension , the operator returns the diagonal matrix
(4) 
For a matrix that is , the operator returns the diagonal matrix
(5) 
The standard Euclidean norm of a vector and Frobenius norm of a matrix are denoted by and , respectively, and defined as follows
(6) 
Two of the applications considered in this paper involve compressible flows governed by the NavierStokes equations. In this case, the fluid state vector is denoted by , and the fluid pressure, density, and velocity magnitude are denoted by , , and , respectively. The freestream conditions are designated by the subscript .
3 Reduction of embedded boundary models
In this work, FSI problems are considered as a vehicle to explain the snapshot compression method proposed for constructing a ReducedOrder Basis (ROB) suitable for the nonlinear reduction of embedded boundary models. It is stressed however that this method is general. In principle, it applies to any EBM framework designed for solving any computational problem. Before presenting this method in Section 4 however, the stage is set here by describing the generic problem of interest, and its reduction process.
3.1 Generic nonlinear fluidstructure interaction problem
Consider the problem of computing an unsteady compressible flow around a rigidly moving body contained in a fixed region . A representative geometry of this problem is illustrated in the left part of Figure 1. Assume that the boundary and fluid/body interface are equipped with the appropriate boundary and fluidstructure transmission conditions, respectively. In an EBM framework, a background Eulerian mesh is typically used to discretize the global domain . The nodes of this mesh that are occluded by at time are usually referred to as ‘‘ghost’’ nodes; the others are referred to as ‘‘real’’ nodes. A semidiscretization on this mesh of the NavierStokes equations governing the fluid subsystem yields a set of nonlinear ordinary differential equations (ODEs) which can be written as
(7) 
where denotes time, , is the number of fluid variables per mesh node, denotes the total number of mesh nodes, and denotes a nonlinear function of that represents the convective and diffusive fluxes.
Without any loss of generality, it is assumed throughout the remainder of this paper that Eq. (7) is discretized in time using an implicit linear multistep scheme. Hence, if denotes a discretization of the timeinterval and , the discrete counterpart of Eq. (7) at timestep is
(8) 
where is the order of accuracy of the chosen timeintegrator and and are two constants that characterize it.
In an EBM framework for CFD, the fluidstructure transmission conditions require special attention because the fluid/body interface does not coincide in general with the nodes of the background mesh. Over the years, a large number of different boundary treatment schemes has been developed for addressing this issue (for example, see Mittal_ARFM_2005 () and Zeng12 () for a review of this and related issues). This is noteworthy because the method proposed in this paper for compressing solution snapshots computed by EBMs based on these schemes is actually independent of them.
3.2 Nonlinear model reduction
In a projectionbased MOR method, the state vector is approximated in a global affine trial subspace as follows
(9) 
where is a matrix whose columns contain the basis of this subspace, is the reduced dimension, denotes the generalized coordinates in this basis at time , and is the initial condition. Substituting Eq. (9) into Eq. (8) yields , which represents an overdetermined system of equations with unknowns. Consequently, the residual is projected onto a test subspace represented by , yielding the square system . In a Galerkin projection, . In a leastsquare projection, Carlberg1 (), where is the Jacobian of (8) with respect to . In either case, the vector of generalized coordinates is then obtained by solving a square system of nonlinear equations using Newton’s method or a preferred variant.
Computing the vector of generalized coordinates using a leastsquares projection approach is equivalent to solving at each timestep the minimization problem Carlberg1 ()
(10) 
For nonlinear, non selfadjoint problems such as those represented in this case by the set of ODEs (7), this approach is more robust than its Galerkin counterpart. Hence, its application in the context of an EBM framework is given here further attention.
Since only a portion of the Eulerian computational domain corresponds to the real flow region, minimizing the norm of the global residual vector  that is, including its ghost component  is neither necessary nor convenient. Instead, let denote a binary vector identifying the fluid/body  or more generally, the ghost/real  partition of the computational domain at time
(11) 
The idea here is that at each timestep, the norm of only the real component of  that is, that associated with the real flow region  needs to be minimized to obtain the vector of generalized coordinates needed to determine the ReducedOrder Model (ROM) approximation (9). In other words, the idea here is that in the context of an EBM, the minimization problem (10) can be reformulated as
(12) 
Equation (12) above is a nonlinear equation. Therefore, it must be solved iteratively using, for example, the GaussNewton method, which can be summarized as
for solve
(13) 
where the superscript designates the th iteration, is the increment of the soughtafter solution at the th GaussNewton iteration, and is determined by a convergence criterion.
Equation (13) constitutes a dimensional GaussNewton ROM of the HighDimensional Model (HDM)  that is, the highfidelity CFD model  represented here by Eq. (7). The normal equations associated with Eq. (13) can be written as
(14) 
where the superscript denotes the transpose.
In any case, whether a Galerkin or leastsquares approach is chosen for determining the vector of generalized coordinates , the timeinvariant ROB must be first determined. In MOR applications, basis functions are usually constructed from solution snapshots. These snapshots are collected, assembled into a matrix, and compressed using, for example, SVD. In an interfacefitted computational framework, this computation is straightforward because the computational domain does not contain a ghost region. In an EBM framework however, this compuation is complicated by the presence and evolution of a ghost/real partition of the computational fluid domain. In the following section, a method for constructing optimal fluid basis functions from snapshots computed by an EBM is presented.
4 Construction of a reducedorder basis for an embedded boundary model
To this effect, this section is organized in three parts. First, the concepts of solution snapshots and a snapshot matrix are recalled. Then, the problem of constructing an optimal ROB for embedded boundary models is formulated as a weighted lowrank matrix approximation problem. Finally, an efficient iterative scheme for solving this problem is specified.
4.1 Eulerian snapshots
A solution snapshot, or simply a snapshot, is defined here as a state vector computed as the solution of Eq. (8) for some instance of its parameters  that is, for some specific time or specific value of the set of flow parameters or boundary/initial conditions underlying this governing equation  on the background Eulerian mesh discretizing the domain . A snapshot matrix is defined as a matrix () whose columns are individual snapshots. For all practical purposes, the main focus of this paper is on unsteady flows and on snapshots associated with different timeinstances . Hence, for .
Because the proposed method for computing a ROB for an embedded boundary model does not utilize information from the occluded region of the computational domain, the ghost components of these snapshots can take arbitrary values. Consequently, the proposed method, which is described below, is independent of the specifics of the EBM framework used to generate the snapshots.
4.2 Weighted lowrank matrix approximation problem
Constructing basis functions from snapshots in the spirit of the POD method  that is, using data compression  can be formulated mathematically as a lowrank matrix approximation problem as follows.
For a given snapshot matrix , find a lower rank matrix that solves the minimization problem
(15) 
where . In this problem, the rank constraint can be taken care of by representing the unknown matrix as , where and , so that problem (15) becomes
(16) 
It is wellknown that the solution of the above lowrank approximation problem is given by the SVD of : specifically, and where .
Unfortunately, when the snapshot matrix is generated using an EBM framework, the solution outlined above cannot be expected to yield an optimal ROB. This is because: (a) the snapshots computed using an EBM contain information from both the flow (or real) region of the computational domain and its occluded (or ghost) region, and (b) this data inconsistency is not accounted for in the standard lowrank matrix approximation problem (15). Hence, this issue is resolved here by proposing the alternative weighted lowrank matrix approximation problem
(17) 
where is a binary mask matrix identifying the ghost/real partition of the snapshot matrix
(18) 
where has already been defined in (11). This binary weighting ensures that the values of at the ghost nodes do not play a role in the construction of the basis functions, which in turn implies that the derived basis functions are optimal for the real components of the snapshots associated with .
In the general case (i.e. ), the weighted lowrank approximation problem formulated above is not reducible to the unweighted problem. Therefore, its solution is not given by the SVD factorization of the snapshots ^{3}^{3}3The weighted lowrank matrix approximation problem is reducible to the unweighted problem only for the special case where ; see B for details.. Furthermore, no closed form solution of this problem is known. Hence, it must be solved by numerical iterations. In this work, the Alternating Least Squares (ALS) algorithm is chosen for this purpose. This algorithm is the most widely applicable and empirically successful approach for solving this and related problems (Chen_2008, ; Okatani_2007, ; Mitra_2010, ).
4.3 Alternating least squares algorithm
The ALS algorithm takes advantage of the bilinearity of the representation . Using the connection between the Frobenius and Eucledian norms and the notation for the Kronecker matrix product (2), problem (17) can be rewritten as
(19a)  
(19b)  
(19c) 
Problem (19b) is a linear leastsquares problem for if is known. Problem (19c) is a linear leastsquares problem for if is known. This suggests the following simple iterative solution procedure:

Guess (for example, initialize as the first leftsingular vectors of )
A more detailed outline of the above algorithm is given in Algorithm 1, where the stopping criterion has been omitted for the sake of brevity (see (Markovsky_book_2012, ; Markovsky_JCAM_2013, ) for a discussion of this topic).
The superscripts of and in Algorithm 1 designate an ALS iteration. The computational cost of this algorithm is roughly times the computational cost of a single thin SVD of the snapshot matrix times the number of iterations for convergence  that is, , where denotes the performed number of iterations. In Section 5, it is shown that in general, suffices to obtain good results. The convergence of this algorithm is monitored using the normalized weighted projection error
(20) 
This error is guaranteed to decrease monotonically to a local minimum (Markovsky_book_2012, ; Markovsky_JCAM_2013, ).
5 Applications
Three FSI problems are considered here to illustrate the method proposed in this paper for constructing basis functions suitable for the nonlinear model reduction of embedded boundary models, and demonstrate its features and performance. The first one is a simple onedimensional model problem based on Burgers’ equation. It has the merit of illustrating the basic idea and being easily reproducible by the interested reader. The second problem focuses on the prediction of a twodimensional unsteady viscous flow around a cylindrical body in largeamplitude heaving motion. Because this problem is formulated in an unbounded fluid domain, it is also solvable using an interface (or body) fitted ALE framework. Hence, it offers a venue not only for assessing the intrinsic performance of the proposed method for constructing a suitable ROB for an embedded boundary model that is more representative than that of the previous example, but also assessing the performance of the resulting nonlinear ROM relative to that of a counterpart ROM constructed using a bodyfitted computational framework. The third problem focuses on the solution of a twodimensional unsteady turbulent flow inside a square cavity containing a rotating ellipsoidal body. In this case, the largeamplitude motion of challenges the robustness if not feasibility of a bodyfitted ALE framework and calls instead for an EBM framework. Hence, this third problem highlights the need for and demonstrates the performance of the computational methodology proposed in this paper.
For each FSI problem outlined above, the governing (fluid) equations are semidiscretized using the central finite difference method. They are discretized in time using the backward Euler implicit scheme with a constant timestep chosen so that does not travel more than one layer of nodes during this timestep. The Newton method is used to solve all nonlinear equations arising from the implicit timediscretization.
Whenever the flow problem of interest is solved using a bodyfitted ALE framework, a ROB for this problem is constructed by computing snapshots of the solution at different time instances , then compressing them using the SVD. In this case, the nonlinear ROM of interest is generated using the GNAT method Carlberg1 (); Carlberg_JCP_2013 () based on a leastsquares projection.
On the other hand, whenever an EBM framework is used to solve the flow problem of interest, a ROB for this problem is constructed by computing snapshots of the solution at different time instances , and compressing them using the ALS algorithm described in Algorithm 1. In this case, the nonlinear ROM of interest is constructed using a variant of the same GNAT method where at each timestep, the vector of generalized coordinates is computed by solving (12).
5.1 Onedimensional fluidstructure model problem based on the viscous Burgers equation
First, the following onedimensional FSI problem based on the viscous Burgers equation, a periodic boundary condition for the fluid subsystem, and nonhomogeneous Dirichlet boundary conditions for the body subsystem in lieu of the fluidstructure transmission conditions is considered
(21) 
where is a Reynoldslike number and is set to , and
(22) 
The above initial boundary value problem models a flow problem in an unbounded, onedimensional, fluid domain in the middle of which a rigid, linear body of length is placed and set into the oscillatory harmonic motion characterized by the amplitude 0.1, frequency , and velocity . At each time , defines the position of the left extremity of , and defines the position of its right extremity. Hence, this problem is a onedimensional instance of the generic FSI problem discussed in Section 3.1.
The simplicity of problem (21) is such that it can be solved using both an ALE framework on a moving mesh, and an EBM framework on a fixed mesh. Hence, both frameworks are considered here, primarily for the purpose of comparisons. In both cases, Eq. (21) is discretized on a uniform Cartesian mesh with ( nodes).
Specifically, due to the unbounded fluid domain assumption and the periodic fluid boundary condition, the bodyfitted ALE framework considered here for the solution of problem (21) is equipped with a rigid motion of the mesh that is identical to that of . As for the considered EBM framework, it is equipped with a firstorder ghost fluid scheme where the ghost values of are populated at each timeinstance only at the two ghost fluid nodes that are the nearest to the left and right boundaries of , using the same value . This EBM framework is illustrated in Figure 2, where the circles filled in black represent the real nodes, the empty ones represent the ghost fluid nodes, and the empty squares designate the subset of ghost fluid nodes whose ghost fluid values are populated. Both computational frameworks deliver the same HDM solution which is graphically depicted in Figure 3.
For each of the bodyfitted ALE and EBM models, three ROBs of dimension , 20, and 40 are constructed using snapshots of the solution of problem (21), and three corresponding nonlinear ROMs are generated.
Figure 4(a) reports on the convergence of the ALS algorithm using the normalized weighted projection error (20). It illustrates a wellknown property of this algorithm, namely, its guarantee to deliver a monotonic convergence to a local minimum. Figure 4(b) reports on the variation with the ROM dimension of the EBM ROM error defined as
(23) 
where is the matrix of snapshot solutions of problem (21), is the ROB computed at the the iteration of the ALS algorithm, and is the matrix of generalized coordinates of the approximate solutions associated with and can be written as
(24) 
 that is, each th column of the matrix product is the EBM ROM solution of problem (21) at time , based on the ROB computed at the th iteration of Algorithm 1.



At the initial (th) ALS iteration, the ROB is constructed by performing the SVD directly on EBM solution snapshots where the ghost values of the solution  that is, the values of this solution at the occluded mesh nodes  are populated with zeros. As shown in Figure 4(b), this approach  referred to here as the ‘‘naive approach’’  delivers a poor performance, which illustrates the need for an alternative method for constructing EBM ROBs. Figure 4(b) also reveals that a single ALS iteration significantly improves the computed EBM ROM solution. The results reported in Figure 4(b) also show that after a single ALS iteration, the constructed EBM ROMs deliver already a comparable accuracy to that of the bodyfitted ALE ROMs of same dimensions.
Figure 5 displays the convergence of the EBM ROM solutions at for , and contrasts these solutions with their counterparts obtained using the bodyfitted ALE ROMs of same dimensions. The reader can observe that both types of ROMs exhibit comparable convergence behavior and accuracy, thereby demonstrating the effectiveness of the proposed method for constructing suitable EBM ROBs.
5.2 Twodimensional fluidstructure interaction problem in an unbounded fluid domain
Next, the simulation of a compressible viscous flow around a heaving rigid cylindrical body with a circular cross section of radius is considered. The fluid is assumed to be a perfect gas with the ratio of specific heats and is modeled using the twodimensional compressible NavierStokes equations. The cylinder is assumed to be infinitely long, and therefore the flow is effectively modeled as a twodimensional problem around a disk. The physical fluid domain is assumed to be unbounded, but the computational fluid domain is bounded by a square of edge length equal to as shown in Figure 6. The freestream Mach number is set to , the Reynolds number based on the cylinder diameter to , and the Prandtl number to .
The center of the disk representing , , is set into the vertical sinusoidal motion
(25) 
where denotes the abscissae of . Hence, the total centertocenter displacement of is , which is a relatively large displacement. Yet, because the fluid domain is also unbounded in this case, and because is rigid, the FSI problem outlined above can be reliably solved by a bodyfitted ALE framework where the CFD mesh is rigidly moved to follow the vertical oscillations of .
Hence, two similar but strictly speaking different meshes are constructed for discretizing the spatial domain shown in Figure 6. The first mesh is designed for the EBM framework. It is nonuniform, has elements whose size in the vicinity of the cylindrical body is uniform and at its finest is characterized by , and embeds the circular body . The second mesh is similar except for the fact that it is a bodyfitted mesh.
The EBM for CFD used for solving this FSI problem is illustrated in Figure 7, where the real fluid nodes are shown as solid circles and the ghost fluid nodes as empty circles. The kinematic transmission condition (or adherence condition) on the fluid/body interface is enforced using a firstorder ghostfluid scheme. At each timestep, the ghost fluid values at the ghost fluid nodes that are the nearest to the fluid/body interface (these are identified by squares in Figure 7), are populated as follows. The ghost fluid velocities at these nodes are set to the velocity of the immersed body . The ghost fluid pressures and densities at these same nodes are computed by constant extrapolation of their counterparts at the nearest neighboring real fluid nodes. On the other hand, the ghost fluid values at all other ghost nodes, which are labelled ‘‘G’’ in Figure 7, are not populated. At the end of each computational timestep, the fluid state vector at the ‘‘transitional’’ nodes (labelled ‘‘T’’ in Figure 7)  that is, the nodes whose status switches from ghost to real  are populated using the ghost fluid values at these nodes from the previous timestep. Snapshots of the vorticity and pressure fluctuations computed using this EBM and the HDM are shown in Figure 8.
First, a pair of unsteady simulations are performed for using the bodyfitted ALE and EBM frameworks outlined above, in order to generate in each case snapshots for model reduction. Then, ROBs and nonlinear ROMs are constructed as outlined above.
Figure 9(a) reports on the convergence of the ALS algorithm using the projection error (20). As expected, this error decreases monotonically.
For this problem, the ROM error is assessed using the relative error in the lift coefficient measured in the norm
(26) 
where and . This error is reported in Figure 9(b). As before, the ROB obtained at the th iteration of the ALS algorithm is computed using the socalled naive approach. Figure 9(b) shows that as expected, the nonlinear ROM constructed using this ROB performs poorly. However, after only one ALS interation, the ROM constructed using the ROB delivers a good performance. The results reported in Figure 9(b) also demonstrate that after a single ALS iteration, the generated EBM ROMs deliver an accuracy that is comparable to that of their ALE ROM counterparts of the same dimensions.
Figure 10 displays the evolution of the instantaneous coefficients of lift and drag, and , where and . The solid curves correspond to computations performed using the EBM HDM, whereas the dashed curves correspond to counterparts performed using the constructed nonlinear ROMs (and a single iteration of the ALS algorithm in the case of the EBM ROMs). The reader can observe that, as expected, the solutions predicted using both sets of ROMs exhibit comparable accuracies and converge to their counterparts obtained using the EBM HDM.



5.3 Twodimensional turbulent fluidstructure interaction problem in a cavity
Finally, the twodimensional computation of a viscous flow inside a square cavity with a rotating ellipsoidal body is considered here. Specifically, the objective is to predict the flow past the ellipsoidal body during the timeinterval in which it rotates a full degrees. This FSI problem with largeamplitude displacements and rotations is chosen because unlike the previous two problems, it cannot be solved using a bodyfitted ALE framework without some form of repeated remeshing to avoid mesh entangling. Hence, it is representative of a large family of FSI problems for which the EBM framework is preferred, if not essential. The geometry of this problem is illustrated in Figure 11.
Again, the fluid is modeled as a perfect gas with the ratio of specific heats . The Reynolds number based on the ellipse tip is set to and the Prandtl number to . The flow is assumed to be initially quiescent.
The square domain is discretized using a uniform Cartesian grid with elements. The same EBM outlined for the previous FSI problem with a heaving cylindrical body is applied to its solution. Computational snapshots of the vorticity contours obtained using the EBM HDM and three different EBM ROMs of dimensions varying between and are illustrated in Figure 12. All EBM ROMs are constructed with a ROB obtained using a single iteration of the ALS algorithm.
As in the previous two cases, the EBM ROM solution is shown to converge to its EBM HDM counterpart when the dimension of the ROM is increased. In particular, the reader can observe in Figure 12(c) that the vorticity field predicted by the EBM ROM of dimension reproduces remarkably well that obtained using the EBM HDM. For , the vorticity fields predicted by the EBM ROM and EBM HDM are almost indistinguishable, thereby demonstrating the effectiveness of the proposed method for constructing ROBs and ROMs for embedded boundary models.
6 Conclusions
Embedded boundary models are popular for the solution of nonlinear problems with evolving domains and/or moving interfaces. The application of projectionbased model order reduction methods to the construction of reducedorder versions of such models calls for the construction first of their underlying reducedorder bases. The collection and compression of solution snapshots using various instances of a computational model of interest is both a popular and effective approach for generating reduced bases. For a traditional interfacefitted computational model, it is a straightforward procedure that has been popularized by the Proper Orthogonal Decomposition and Singular Value Decomposition methods. For an embedded boundary computational framework, this approach faces the problem of dealing with the partitioning of the overall computational domain into a real subdomain and a ghost one corresponding to the region of the overall computational domain that is occluded by the embedded boundary. For timedependent problems, this complication is further exacerbated by the evolution in time of this partitioning as it implies an evolution in time of the partitioning between meaningful (real) and meaningless (ghost) entries of the collected solution snapshots. This makes the problem of compressing these snapshots in view of constructing a reducedorder basis a difficult task. To address this issue, this paper formulates the snapshot compression problem as a weighted lowrank approximation problem where the binary weighting identifies the evolving component of the individual simulation snapshots, and proposes to solve this problem using the Alternating Least Squares (ALS) algorithm. This approach is applicable in principle to any embedded boundary value problem. It is successfully demonstrated in this paper for three different fluidstructure interaction problems of increasing complexity. In all considered cases, it is shown to deliver reducedorder bases and models that effectively reproduce the highdimensional solutions even when the flow is vortexdominated, and the immersed body underoges large displacements and rotations.
Appendix A MATLAB implementation of the ALS algorithm
The following is a simple MATLAB implementation of Algorithm 1 presented in Section 4.3.
The inputs of the function ALS
are the snapshot matrix X
where the ghost values are set to zero, a binary
weighting matrix M
, the desired dimension of the subspace of approximation k
, and the maximum number of
iterations it_max
. The outputs of this function are the lowrank approximation matrices U
and V
.
Appendix B Solution of the weighted lowrank matrix approximation problem when
If is a rank 1 weighting matrix, it can be written as for some and . In this case,
(27a) 
Defining as , as , and as , the above result can be rewritten as
(28a) 
where . This shows that when is a rank 1 weighting matrix, the weighted lowrank approximation problem (17) reduces to its unweighted counterpart (15).
An example of a rank 1 weighting matrix is the weighting matrix appropriate for an embedded boundary model where the embedded body is stationary  that is,  but the flow is unsteady. In this case,
(29) 
Acknowledgments
The authors acknowledge partial support by the Office of Naval Research under Grant N000141110707, and partial support by the Army Research Laboratory through the Army High Performance Computing Research Center under Cooperative Agreement W911NF0720027. The content of this publication does not necessarily reflect the position or policy of any of the aforementioned institutions, and no official endorsement should be inferred.
References
References
 (1) J. Donea, An arbitrary LagrangianEulerian finite element method for transient fluidstructure interactions, Computer Methods in Applied Mechanics and Engineering 33 (1982) 689723.
 (2) C. Farhat, M. Lesoinne, N. Maman, Mixed explicit/implicit time integration of coupled aeroelastic problems: threefield formulation, geometric conservation and distributed solution, International Journal for Numerical Methods in Fluids 21 (1995) 807835.
 (3) C. Farhat, A. Rallu, K. Wang, T. Belytschko, Robust and provably secondorder explicitexplicit and implicitexplicit staggered timeintegrators for highly nonlinear fluidstructure interaction problems, International Journal for Numerical Methods in Engineering 84 (2010) 73107.
 (4) R. Loubère, P.H. Maire, M. Shashkov, J. Breil, S. Galéra, ReALE: A reconnectionbased arbitraryLagrangianEulerian method, Journal of Computational Physics 229 (2010) 47244761.
 (5) C. Farhat, C. Degand, B. Koobus, M. Lesoinne, Torsional springs for twodimensional dynamic unstructured fluid meshes. Computer Methods in Applied Mechanics and Engineering 163 (1998) 231245.
 (6) C. Degand, C. Farhat, A threedimensional torsional spring analogy method for unstructured dynamic meshes. Computers and Structures 80 (2002) 305316.
 (7) M. Fossati, R.A. Khurram, W.G. Habashi, An ALE mesh movement scheme for longterm inflight ice accretion, International Journal for Numerical Methods in Fluids 68 (2012) 958976.
 (8) C. S. Peskin, The immersed boundary method, Acta Numerica 11 (0) (2002) 479517.
 (9) R. Mittal, G. Iaccarino, Immersed boundary methods, Annual Review of Fluid Mechanics 37 (2005) 239261.
 (10) K. Wang, A. Rallu, J. F. Gerbeau, C. Farhat, Algorithms for interface treatment and load computation in embedded boundary methods for fluid and fluidstructure interaction problems, International Journal for Numerical Methods in Fluids 67 (9) (2011) 11751206.
 (11) C. Farhat, J.F. Gerbeau, A. Rallu, FIVER: a finite volume method based on exact twophase riemann problems and sparse grids for multimaterial flows with large density jumps, Journal of Computational Physics, 231 (2012) 63606379.
 (12) M. Rewieński, J. White, Model order reduction for nonlinear dynamical systems based on trajectory piecewiselinear approximations, Linear algebra and its applications 415 (2) (2006) 426454.
 (13) T. Lieu, C. Farhat, M. Lesoinne, Reducedorder fluid/structure modeling of a complete aircraft configuration, Computer Methods in Applied Mechanics and Engineering 195 (2006) 57305742.
 (14) D. Amsallem, C. Farhat, An interpolation method for adapting reducedorder models and application to aeroelasticity AIAA Journal, 46 (2008) 18031813.
 (15) P. Astrid, S. Weiland, K. Willcox, T. Backx, Missing point estimation in models described by proper orthogonal decomposition, Automatic Control, IEEE Transactions on 53 (10) (2008) 22372251.
 (16) S. Chaturantabut, D. Sorensen, Nonlinear model reduction via discrete empirical interpolation, SIAM Journal on Scientific Computing 32 (5) (2010) 27372764.
 (17) K. Carlberg, C. BouMosleh and C. Farhat, International Journal for Numerical Methods in Engineering, 86 (2) (2011) 155181.
 (18) M. J. Balajewicz, E. H. Dowell, B. R. Noack, Lowdimensional modelling of highreynoldsnumber shear flows incorporating constraints from the NavierStokes equation, Journal of Fluid Mechanics 729 (2013) 285308.
 (19) K. Carlberg, C. Farhat, J. Cortial, D. Amsallem, The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows, Journal of Computational Physics 242 (2013) 623647.
 (20) J. Du, F. Fang, C.C. Pain, I.M. Navon, J. Zhu, D.A. Ham, POD reducedorder unstructured mesh modeling applied to 2D and 3D fluid flow, Computers and Mathematics with Applications, 65 (2013) 362379.
 (21) D. Xiao, F. Fang, A. G. Buchan, C.C. Pain, I.M. Navon, J. Du, G. Hu Nonlinear model reduction for the NavierStokes equations using residual DEIM method, Journal of Computational Physics, 10.1016/j.jcp.2014.01.011.
 (22) N. Srebro, T. Jaakkola, Weighted lowrank approximations, in: ICML 3 (2003) 720727.
 (23) A. M. Buchanan, A. W. Fitzgibbon, Damped newton algorithms for matrix factorization with missing data, in: Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, IEEE 2 (2005) 316322.
 (24) P. Chen, Optimization algorithms on subspaces: Revisiting missing data problem in lowrank matrix, International Journal of Computer Vision 80 (1) (2008) 125142.
 (25) D. L. Donoho, Compressed sensing, Information Theory, IEEE Transactions on 52 (4) (2006) 12891306.
 (26) E. J. Candès, The restricted isometry property and its implications for compressed sensing, Comptes Rendus Mathematique 346 (9) (2008) 589592.
 (27) E. J. Candès, B. Recht, Exact matrix completion via convex optimization, Foundations of Computational Mathematics 9 (6) (2009) 717772.
 (28) J. F. Cai, E. J. Candès, Z. Shen, A singular value thresholding algorithm for matrix completion, SIAM Journal on Optimization 20 (4) (2010) 19561982.
 (29) P. Jain, P. Netrapalli, S. Sanghavi, Lowrank matrix completion using alternating minimization, in: Proceedings of the 45th annual ACM symposium on Symposium on theory of computing, ACM (2013) 665674.
 (30) B. Vandereycken, Lowrank matrix completion by riemannian optimization, SIAM Journal on Optimization 23 (2) (2013) 12141236.
 (31) T. Okatani, K. Deguchi, On the wiberg algorithm for matrix factorization in the presence of missing components, International Journal of Computer Vision 72 (3) (2007) 329337.
 (32) K. Mitra, S. Sheorey, R. Chellappa, Largescale matrix factorization with missing data under additional constraints, in: Advances in Neural Information Processing Systems (2010) 16511659.
 (33) I. Markovsky, Low rank approximation: algorithms, implementation, applications, communications and control engineering, Springer (2012).
 (34) I. Markovsky, K. Usevich, Software for weighted structured lowrank approximation, Journal of Computational and Applied Mathematics 256 (2013) 278292.
 (35) X. Zeng, C. Farhat, A systematic approach for constructing higherorder immersed boundary and ghost fluid methods for fluid and fluidstructure interaction problems, Journal of Computational Physics 231 (2012) 28922923.