key: cord-0830348-gbqsvarv
authors: Viguerie, Alex; Barros, Gabriel F.; Grave, Mal'u; Reali, Alessandro; Coutinho, Alvaro L.G.A.
title: Coupled and Uncoupled Dynamic Mode Decomposition in Multi-Compartmental Systems with Applications to Epidemiological and Additive Manufacturing Problems
date: 2021-10-12
journal: Computer Methods in Applied Mechanics and Engineering
DOI: 10.1016/j.cma.2022.114600
sha: 1c30bc70b627806667ac35a70a8f5350e8c2634e
doc_id: 830348
cord_uid: gbqsvarv

Dynamic Mode Decomposition (DMD) is an unsupervised machine learning method that has attracted considerable attention in recent years owing to its equation-free structure, ability to easily identify coherent spatio-temporal structures in data, and effectiveness in providing reasonably accurate predictions for certain problems. Despite these successes, the application of DMD to certain problems featuring highly nonlinear transient dynamics remains challenging. In such cases, DMD may not only fail to provide acceptable predictions but may indeed fail to recreate the data in which it was trained, restricting its application to diagnostic purposes. For many problems in the biological and physical sciences, the structure of the system obeys a compartmental framework, in which the transfer of mass within the system moves within states. In these cases, the behavior of the system may not be accurately recreated by applying DMD to a single quantity within the system, as proper knowledge of the system dynamics, even for a single compartment, requires that the behavior of other compartments is taken into account in the DMD process. In this work, we demonstrate, theoretically and numerically, that, when performing DMD on a fully coupled PDE system with compartmental structure, one may recover useful predictive behavior, even when DMD performs poorly when acting compartment-wise. We also establish that important physical quantities, as mass conservation, are maintained in the coupled-DMD extrapolation. The mathematical and numerical analysis suggests that DMD may be a powerful tool when applied to this common class of problems. In particular, we show interesting numerical applications to a continuous delayed-SIRD model for Covid-19, and to a problem from additive manufacturing considering a nonlinear temperature field and the resulting change of material phase from powder, liquid, and solid states.

In recent years, data-driven techniques have become ubiquitous in science and engineering, as large quantities of high-dimensional data become more available and easier to handle on modern computer architectures [13] . Though algorithms and techniques for the numerical simulation of physical systems have improved, such simulations still demand considerable resources in terms of both time and computing power. For such problems, the emerging field of Scientific Machine Learning (SciML) represents a potential way to reduce such costs by exploiting the inherent structure of such large datasets. Such techniques can be used in lieu of simulation or as a tool to extract useful information from numerical simulation data. SciML can be useful, for instance, to accelerate convergence, extract long-term behavior using only short-term simulations, or reduce the number of forward solves in applications involving the solution of hundreds, or possibly thousands, of numerical simulations (such as those common in optimization or uncertainty quantification).

There is a wide range of SciML methods, which have different theoretical bases and intended use cases [14] . For more information on the field of SciML as a whole, we kindly refer the reader to [51, 24, 45, 14, 39, 40] . In the present work, we discuss Dynamic Mode Decomposition (DMD), an unsupervised SciML method that exploits the underlying low-rank structure in many large datasets in order to extract only the most relevant information. This low-rank reconstruction may then be employed for diagnostic purposes or to make future state predictions. While DMD may be used on both numerical simulations [49] and experimental data [50] , in the present work, we restrict our attention to the former, though the discussion applies equally well to experimental data.

The standard DMD procedure takes m snapshots of state measurements, and uses them to build a linear operator that best approximates the nonlinear system dynamics. While this idea has a strong theoretical motivation, in practice, DMD often fails to provide adequate results for extrapolation and prediction. This DMD behavior is particularly common for problems exhibiting a highly nonlinear transient behavior. In such cases, while DMD may still provide valuable diagnostic information [13, 3, 44, 32] , it may fail to provide an acceptable predictive behavior, in some cases failing even to provide an accurate reconstruction of the training period. Attempts to overcome these shortcomings have been addressed in, e.g., [48, 17, 41] , in which various different methods, including statistical procedures, time delay embeddings, and augmenting the observable space with additional nonlinear variables have been proposed to improve such performance.

Types of system for which DMD may fail to provide acceptable predictions are those of reaction-diffusion type, commonly found in biological and physical sciences [41, 4, 55] . Such systems often feature a compartmental structure, in which the total mass of the system is decomposed into various compartments (or species) [54] . In the special case of closed compartments, the total mass of the system remains constant, with no additional sources or sinks. For large scale-systems, DMD may be performed compartment-wise, applying the technique on each physical compartment separately [4] . Such an approach is understandable in the presence of many compartments and large-scale numerical discretizations, owing to the reduced costs of performing such discretizations separately, possibly in parallel. However, while this method may still provide acceptable results for diagnostics and feature identification, it can provide suboptimal results when used for prediction, or even simply reconstruction, of the desired system. In particular, many important system-wide characteristics, such as conservation of physical quantities across compartments, may no longer be respected [4] .

In the present work, we advocate a fully-coupled approach; that is, DMD should be performed on the entire compartmental system, rather than on the single compartments individually. Applying the technique in this manner may provide an advantage in system reconstruction and prediction, as knowledge of even a single compartment is inherently linked to the evolution of the other compartments. When performing DMD in this manner, the underlying structure governing the entire system may be properly taken into account and reconstructed, even when extrapolating the system in time to predict future states. Notably, for problems featuring a closed-compartment structure, in which a total quantity of the system remains constant in time, this property is maintained in the DMD reconstruction. Thus, for such a class of problems, we show that DMD, when applied properly, is particularly well-suited to provide reasonable forecasts. This paper is structured as follows: Section 2 describes the relation between the discretization of multi-compartment reaction-diffusion PDEs in space and time and the DMD algorithm. Section 3 describes the proposed methodology of performing DMD on the fully coupled system, rather than on the individual compartments, and provides some theoretical results justifying this choice. In Section 4 we describe the numerical applications in this study: the use of DMD on a continuous delayed-SIRD model from epidemiology, and a model from additive manufacturing which considers a nonlinear temperature field and the resulting change of material phase from powder, liquid, and solid states. We demonstrate how compartment-by-compartment reconstructions both fail to provide useful predictions and violate the conservation of mass. We then show that performing DMD on the fully coupled system eliminates these problems, allowing for the recovery of useful extrapolations while respecting the problem physics. In Section 5, we draw our final remarks and conclusions.

The use of numerical methods for the approximation of partial differential equations (PDEs) is widely considered in many research topics and industrial applications. Many natural and artificial phenomena can be mathematically modeled using PDEs, and, by solving these equations, we can describe and analyze multiple processes that help improving the society as a whole, from finances and social behavior to complex physical and engineering systems. Many numerical methods such as , e.g., the finite difference method, the finite element method, and the boundary element method are reliable mathematical tools for the approximation of PDEs, allowing a proper quantification of the different quantities of interest.

In the present study, we focus on using the finite element method for the spatial discretization of a system of compartmental reaction-diffusion equations. That is, consider a system of generic transient parametric PDEs for a vector u = [u 1 , u 2 , ..., u n ] T , where the u i correspond to the vector compartments. We define ζ as a vector of parameters governing the control of the nonlinear system, and f = [f 1 (t; ζ), f 2 (t; ζ), ..., f n (t; ζ)] T as given functions. We then consider the equation:

where the subscripts i, j indicate that the operator N i,j defines the interaction between the i-th and j-th compartments. The equation is equipped with boundary and initial conditions for each compartment

where g = [g 1 , g 2 , ..., g n ] T and h = [h 1 , h 2 , ..., h n ] T are the vectors containing the Dirichlet and Neumann boundary conditions for each equation, respectively, T is the final time, and u 0 == [u 01 , u 02 , ..., u 0n ] T are the initial conditions for u. The equations are solved on a bounded domainΩ composed by the domain Ω ⊂ R nsd and Lipschitz continuous boundaries Γ D ∪ Γ N = Γ ⊂ R nsd−1 . Also, n is the unit outward normal to Γ N . Considering the classical finite element method for the spatial discretization of the PDEs, the bounded domain is split into a finite number of elements with non-overlapping domain Ω e ⊂ R nsd and boundaries Γ e ⊂ R nsd−1 . Moreover, the finite element method requires the equations to be converted in their weak forms by integrating Eq. (1) against a weighting function w ∈ H 1 (Ω) and, then, applying the divergence theorem on the resulting equation. The weighting functions exist on H 1 (Ω) which is the Sobolev space of the square-integrable functions with an integrable first weak derivative, allowing the weak form to be continuous. Considering P k (Ω e ) the space of polynomials of degree equal or less than k over Ω e , the function spaces are defined as:

That said, the weak formulation for the spatial discretization for the nonlinear system described on Eq. (1) is: Find

where the L 2 inner product over the domain Ω is indicated by (·, ·). The choice for number of elements, element type, finite element formulation (i.e., Galerkin, stabilized formulations, etc.) is of major importance and depends largely on the nature of the problem to be solved. For instance, if the reaction term on the diffusion-reaction equations is dominant, the standard Galerkin formulation may fail to produce physical results, being well replaced by stabilized formulations such as the SUPG [12] or Variational Multiscale methods [27, 46, 1, 16, 5] . For simplification regarding the notation, the superscript "h" will be omitted as the discretization of continuous variables has been applied.

Since the system of equations has been partially approximated, temporal discretization of Eq. (5) is required for the complete resolution of the problem. From this point, the equation can be translated into a discrete-time dynamical system. In this system, the state vector u at the time instant k + 1 can be written such that:

that is, the system can now be considered as a temporal update of the state vector u k+1 containing the information of all the compartments on a given instant at a discrete level. The matrix F is the discrete-time flow map of the system and contains information regarding the parameters ζ, mesh size, solver properties, etc. The construction of u k+1 is done sequentially with the information of the previous state vector. Considering the evolution in time of a discretized PDE as a dynamical system is fundamental for the DMD application on numerical simulations.

DMD is a data-driven method able to extract the most dynamically relevant structures (DMD modes) on a spatiotemporal dataset, where each DMD mode is associated with frequency, and amplification/damping terms [36] . The coherent structures are related in terms of space and time. The method is completely equation-free and data-driven, meaning that little to no assumptions on data must be considered for its applications. DMD was first applied in fluid dynamics applications [47, 49] , being expanded to many other applications such as epidemiology [44, 4] , urban mobility [2] , biomechanics [15] , climate [37] and aeroelasticity [20] , especially in structure extraction from data and control-oriented methods. The method consists of creating a linear map of the dynamics of a given spatio-temporal dataset, even if the dynamics is nonlinear, by projecting the finite-dimensional nonlinear data using an infinite dimension operator able to linearly represent the flow map present on (6) for all time steps.

Looking at transient PDEs through the lens of dynamical systems is crucial for the application of DMD and the system described on Eq. (6) can be considered for this purpose. Let Y j be a dataset containing the state vectors for the compartment j resulting of the dynamical system u j for k = 0, 1, . . . , m, where m + 1 is the total number of entries. The dataset

can be rearranged into two matrices

where n is the number of nodes in the mesh. DMD aims at constructing the best fit approximation of the linear operator A that maps dataset Y into dataset Y , that is,

The matrix A can be computed as A = Y Y † , where Y † is the Moore-Penrose pseudoinverse of Y . However, we avoid the computation of the full matrix A since it is a n × n matrix. Also, the computation of the full Moore-Penrose pseudoinverse is not advisable due to its ill-conditioning. Instead, we can compute the SVD of Y as

where U ∈ R n×m and V ∈ R m×m are the left and right singular vectors and Σ ∈ R m×m is a diagonal matrix with real, non-negative, and decreasing entries named singular values. The singular values σ 0 ≥ σ 1 ≥ σ 2 ≥ · · · ≥ σ m−1 are hierarchical and can be interpreted in terms of how much the singular vectors influence the original matrix Y . For the DMD procedure, considering the Eckart-Young Theorem [18] , the optimal low-rank update approximation matrix Y , when subjected to a truncation rank r, can be written as

where U r ∈ R n×r is a matrix containing the first r columns of U, V r ∈ R m×r contains the first r columns of V, and Σ r ∈ R r×r is the diagonal matrix containing the first r singular values. The pseudoinverse can be approximated as

and, instead of computing A ∈ R n×n , we can obtainÃ, a r × r projection of A, as

Note thatà is unitarily similar to A. Further mathematical details regarding the optimization problem (the best-fitting matrixÃ) and the influence of the Eckart-Young Theorem on constraints of the problem can be found in [28] . Now we can compute the eigendecomposition ofÃ:à W = WΛ, (13) where Λ is a diagonal matrix containing the discrete eigenvalues λ j and the matrix W contains the eigenvectors φ j of A. The DMD basis can be written as:

and the discrete signal reconstruction as:

The continuous counterpart to (15) reads:

with b being the vector containing the projected initial conditions such that b = Ψ † u 0 , and Ω eig is a diagonal matrix whose entries are the continuous eigenvalues ω i = ln(λ i )/∆t o , where ∆t o is the time step size between the outputs. The form of (16) can be regarded as a generalization of the Sturm-Liouville expansion for a differential problem:

where ψ i and ω i are the i-th Sturm-Liouville eigenfunctions and eigenvalues for a given differential operator.

Regarding the use of DMD on systems of the type (1), one may employ one of two approaches. The first approach, which we will refer to as the uncoupled approach, involves treating each u i separately in the DMD process. That is, we define the snapshot matrices:

We then perform the DMD algorithm n times, obtaining n separate reconstructions, one for each compartment. This approach was used in [4] , where it was shown that it does not always yield acceptable results. In particular, in some situations, we may observe large spurious oscillations and rapid deterioration of error [4] , and, in extreme cases, it may even fail to reproduce the dataset used for the training, as has been observed in many DMD applications [48] .

In contrast, we may consider the following alternative:

The DMD reconstruction is then performed once, on a larger system. To recover the signal reconstruction of the i-th compartment for a system containing n compartments of dimension s, we recall that the DMD reconstruction has the same topology as the input data, and hence:

However, it is also the case that, even for large-scale PDE problems, the compartments u i are not, in fact, independent quantities; rather, they may interact with one another in a complex, nonlinear fashion. With this in mind, it makes sense that a proper reconstruction of a given compartment may require adequate knowledge of all of the other compartments in a nonlinear system.

In terms of computational aspects, the uncoupled or coupled approach should be considered carefully. While the overall amount of data to be processed is identical for both methods, computational and memory costs are not. For instance, loading the complete dataset when applying the coupled DMD for a sufficiently large number of compartments could be prohibitive on a given system, depending on the number of degrees of freedom and available memory. However, if this is the case, the uncoupled approach may alleviate these restrictions by loading one (or a batch of) compartment at a time from disk. As for computational time, the algorithmic complexity of DMD is primarily driven by the SVD algorithm. This dependence is inherited when evaluating the use of the uncoupled or coupled cases. The classical SVD algorithm [26] computes the singular values of a dense matrix at the cost of O(n 3 ), implying that the coupled approach will scale cubically with the number of compartments, while the uncoupled approach will scale linearly. This bad scaling leads to a prohibitively high cost of the coupled approach for a large n and/or a large number of compartments.

However, it is worth noting that the classical SVD algorithm considers the complete factorization of the dataset, while DMD only relies on the first k singular values and vectors. In a situation with a large n or a large number of compartments, the extra cost of the coupled DMD may be circumvented by the Randomized Singular Value Decomposition (rSVD) [25, 19] , which computes the approximate rank-k SVD with a computational complexity of O(nmk). Applying the rSVD may enable coupled/uncoupled approaches to achieve a similar overall complexity. We note that there are many additional computational aspects to consider, but such an in-depth discussion is beyond the scope of the present work. We also note that the uncoupled algorithm lends itself more naturally to parallelization. We do acknowledge, however, that such aspects are important and certainly worthy of future investigation.

We now present three mathematical results suggesting that the coupled approach may provide more accurate and physically consistent reconstructions than the uncoupled one for closed-compartment systems. The first two theorems can be expected to hold for all such systems studied in the present work. The third requires more restrictive conditions which, while not guaranteed to be satisfied, we have found in practice to be often satisfied enough for the result to remain useful. Theorem 1: For a compartmental dynamical system u = [u 1 , u 2 , ..., u n ] T such that:

the exact coupled DMD reconstruction u = [ u 1 , u 2 , ..., u n ] T is also such that:

Proof. From (21), we may conclude that, for each instance k, and denoting the column vector of all ones as 1,

for some constant C 1 . This implies that, for (19) ,

and in particular,

We then have that:

This establishes the base case, for which one may apply mathematical induction for u k using the exact argument above. We note that, in practice, we do not expect this condition to hold exactly, as we use approximations to A. However, we should observe this behavior, at least approximately, in the coupled reconstruction. Although we obtained the result for the discrete time-flow map, we also expect it to hold for general t. A proof of this is desired.

Theorem 2. For problems in which some subset of compartments satisfies the hypotheses of Theorem 1, Theorem 1 holds on the DMD reconstruction of the subset.

Proof. As the ordering of compartments is arbitrary, reorder Y such that compartments 1, , 2, ..., j are the first j compartments. Then repeating the arguments of Theorem 1 with

instead of 1, establishes that the closed-compartment property holds on the subset.

It remains to be shown that the compartment reordering only results in a similar reordering on the DMD reconstruction, and does not affect the reconstruction dynamics. Let Y denote the snapshot matrix in the original ordering of the system and Y the reordered snapshot matrix. Note that: Y = PY, where P is a permutation matrix. Then,

Hence, A and A are unitarily similar, and have the same eigenvalues. Further, for a given eigenvector w of A, Pw is an eigenvector of A. The DMD modes are therefore identical up to the reordering permutation, yielding the desired result.

Theorem 3. For a compartmental problem which satisfies the hypotheses of Theorem 1, if we have additionally that u i ≥ 0, u i ≥ 0 at all (x, t) for each i, then the L 1 error of the DMD reconstruction:

Proof. From the assumption of non-negativity on the compartments, it follows that:

We then observe that:

Remark 1. Systems which satisfy the hypotheses of Theorems 1 and 2 are quite common, making these theorems valuable tools. Unfortunately, the non-negativity hypothesis on the u i for Theorem 3 does not hold in general, and these quantities may become negative. In practice, we have found that this assumption remains satisfied for many problems, at least for a moderate period of time after the end of the training period. In such a regime, we may expect the error bound given by Theorem 3 to hold, preventing error blowup. These results suggest, however, that a variant of DMD guaranteeing non-negativity of u i could be a significant development and is a direction worthy of investigation, as it would establish DMD as a valuable tool for predicting the dynamics of closed compartmental systems.

Remark 2. It is common in DMD to augment the variable space with nonlinear observables, typically functions of the state variables, to better capture the nonlinear dynamics of the system. Such an augmentation is motivated by the connection between DMD and the Koopman operator. We note that the result of Theorem 2 holds equally well if the variables which do not obey the closed-compartment structure are not necessarily compartments themselves but additional nonlinear variables. More concretely, one may consider a two-compartment model in which

at each time-step. If one wishes to perform DMD on an augmented variable space [u 1 , u 2 , g(u 1 , u 2 )], with g some nonlinear function depending on u 1 and u 2 , we would still expect that

despite the presence of the nonlinear variable.

In this section, we seek to understand the differences in accuracy between the coupled and uncoupled DMD formulations for both reconstruction as well as extrapolation and prediction. We select two different compartmental models in order to assess the effectiveness of the two approaches:

1. A delayed-PDE SIRD model for COVID-19, which separates the total population into compartments based on their infection and immunity status; 2. A problem inspired by additive manufacturing, consisting of a nonlinear thermal PDE coupled with a three-ODE compartmental model consisting of three differential equations in space, corresponding to the powder, liquid, and solid phases.

We have selected these examples due to their different physics, application areas, and dynamics. They share, however, a mass-conservative closed-compartmental structure. Further details, including the formal mathematical definitions of each model, are given in the following.

The outbreak of COVID-19 in 2020 has led to a surge in interest in the mathematical modeling of infectious diseases. Disease transmission may be modeled by compartmental models, in which the population under study is divided into compartments and has assumptions about the nature and time rate of transfer from one compartment to another [10] , owing to the original SIR framework of Kermack and McKendrick [31] 2 . While less common than the ordinary differential equation (ODE) version of such models, compartmental models for epidemic spread spatial transmission through the use of PDE models have also been applied to study both the COVID-19 epidemic as well as other diseases [54, 53, 23, 22, 8, 7, 9, 29, 30, 43] .

Here, we work with a spatio-temporal delayed SIRD model, related to that presented in [23] and given by the following set of equations:

where s(x, t), i(x, t), r(x, t), and d(x, t) denote the densities of the susceptible, infected, recovered, and deceased populations, respectively. The sum of all compartments with the exception of d(x, t) is represented by n pop which is the total living population. β i and β e denote the transmission rates between symptomatic and susceptible individuals and presymptomatic and susceptible individuals, respectively (units days −1 ), γ the symptomatic recovery rate (units days −1 ), δ represents the mortality rate (units days −1 ), and ν s , ν i , ν r are the diffusion parameter matrices of the different population groups as denoted by the sub-scripted letters (units km 2 persons −1 days −1 ). Note that all these parameters can be considered time and space-dependent. The time-delay σ incorporates the lag effect present between the time of infection and development of symptoms, and seeks to model the delay observed in practice when new restrictions or other measures are introduced.

we may write the system in the generic form given by (1) . For the numerical solution of (30)-(33), we discretize in space using a Galerkin finite element variational formulation. The resulting systems of equations are stiff, leading us to employ implicit methods for time integration. We apply the Backward Euler formula (BE), which offers first-order accuracy while remaining unconditionally stable. This choice is motivated by the ease of incorporating this time-stepping scheme with the presence of the delay term, which is known to make the use of multi-step methods more complex [6] . We implement the whole model in libMesh [33] . One may find more details about the methods in [53, 21] .

We now use the model (30)- (33) to simulate the outbreak of COVID-19 in the U.S. state of Georgia. The variables ν and β i,e vary in time, and are shown in Fig. 1 [22] . The other parameters are assumed fixed for simplicity, although an even more realistic model may incorporate time-dependence for these parameters as well.

We let δ = 1/180 days −1 , γ=1/24 days −1 and σ = 7 days. It was shown in [23] that under the condition:

we expect asymptotically stable behavior from the model (30)- (33) , and further under the more restrictive condition

the compartments of the equations (30)- (33) should remain positive, preventing nonphysical oscillations into negative values. A straightforward calculation shows that, for the adopted values of γ, δ, and σ, we may expect both stable and physical behavior.

The physical domain, mesh, and population distributions were reconstructed following the procedure outlined in [22] . As this is a delay-differential equation model, we let i(t) = i 0 for t ∈ [−σ, 0]. We use a fixed time-step of ∆t = 0.25 days. Since σ = 7 days is a fixed, exact multiple of the time-step size, we may directly use previously-computed solutions of i for i(t − σ), and do not need to use more sophisticated interpolation techniques. We note that extensions of (30)-(33) incorporating adaptive time-stepping and/or a time-and state-dependent σ will require a more involved computational procedure for i(t − σ). We considered an unstructured triangular mesh consisting of 22310 elements, and used a piecewise linear finite element approximation. In Fig. 2 , we plot the results of the measured fatalities in Georgia against the simulated data. We observe excellent agreement overall, with an R 2 of 99.52% amd a relative L 2 norm error of 5.74%. We note that there is a jump in the measured data; this is due to previous unaccounted-for fatalities being retroactively added to the dataset. The evolution of the infected compartment i over several different phases of the epidemic is shown in Fig. 3 .

Having validated the baseline simulation at a satisfactory level, we now analyze the DMD reconstructions. We perform both the uncoupled and coupled DMD. In both cases, we train the DMD reconstruction using the first 250 days of the simulation, and then extrapolate over the following 25 days. Thus, the window for extrapolation is 10% of the training period, consistent with the analysis performed in [4] .

We first confirm Theorem 1 numerically. As we do not consider new births or non-COVID deaths in (30)-(33), we therefore expect that:

for all t, representing a conservation of total mass in the system. We compute the quantity Ω (s + i + r + d)dΩ at all t for both the uncoupled and coupled reconstructions and plot the results in Fig. 4 .

As we can see in the plot, the coupled DMD reconstruction is mass-conservative, for both the reconstruction and extrapolation periods, confirming the theoretical result. In contrast, the uncoupled DMD reconstruction fails to conserve mass, even failing during the training period. This provides a strong justification for the use of the coupled reconstruction formulation, as it is able to properly represent the underlying physics of the model.

We now turn our attention to the error behavior, defined as the relative L 2 error between the DMD reconstructions and the computed compartment for the reference simulations. The results are shown in Fig. 5 . As the figure shows, the uncoupled simulation is unable to provide meaningful results for the infected, deceased, and recovered compartments, performing only an acceptable reconstruction for the susceptible compartment. In contrast, the coupled reconstruction gives a good reconstruction of each compartment over the training period, with errors consistently lower than 10%, and Figure 4 : Comparison in conservation of mass between the coupled and uncoupled DMD reconstructions. We observe that the coupled reconstruction successfully conserves mass, even past the training period, while the uncoupled reconstruction does not respect this important property. Figure 5 : Error behavior in time for Georgia simulation. We see that the uncoupled DMD reconstructions fail for several compartments, giving unacceptable results. In contrast, the coupled DMD formulation is able to properly reconstruct the training dynamics, as well as provide reasonable extrapolations over a three-week interval, with the error compared to the reference simulation less than 1% in several compartments as far as 15 days from the end of the training period.

also for a significant projection window. For the infected compartment, the error between reconstruction and reference is less than 1% as far as 15 days from the end of the training period (Fig. 6) . This is consistent with the theory, as the mass-conservation suggests that, as long as our reconstructions remain non-negative (as is the case here), the error in each compartment will remain under control.

Overall, we observe that the coupled DMD reconstruction provides reasonably accurate predictions over a short time interval, and is able to respect the underlying physics of the compartment. Indeed, even as the solution quality begins to deteriorate as we move further from the training period, the mass of the system remains conserved, suggesting that Figure 6 : The infected compartment for Georgia, with the reference simulation (left) compared with the DMD reconstruction (right), 15 days after the end of the training period. We observe excellent agreement between the two quantities, with the relative L 2 error between them less than one percent.

the overall system physics is well-incorporated into the DMD reconstruction. On the other hand, the compartmentby-compartment uncoupled DMD reconstruction neither respects the mass-conservation nor provides acceptable reconstructions and/or predictions. For this example, we see that the coupled DMD dynamics may indeed be applied for both signal reconstruction and short-term predictions, even when uncoupled DMD performs poorly.

In this section, we consider a two-dimensional problem inspired by additive manufacturing which comprises a layer of powder heated by a laser. As the powder reaches its melting point, it becomes a liquid, whereby upon cooling it becomes a solid. When reheated, the solid may again become liquid; however, the change from powder to liquid is irreversible.

Mathematically, we consider the model used in [52] , which consists of a partial differential equation for the temperature T in a domain Ω and three ordinary differential equations at each point in Ω for the quantities φ p , φ l , and φ s : the phase fractions of the powder, liquid, and solid quantities respectively. The equations then read:

complemented by suitable initial and boundary conditions. T s is the temperature where melting begins, and T f the temperature where melting is complete. ρ is the material density, c p the specific heat, and κ the thermal conductivity. These parameters depend on both the temperature field and material phase. This dependence is understood if not explicitly denoted. For these values, we base our simulation on stainless steel 316L, a common alloy used in additive manufacturing. More information about the material parameters and their specific relationship to temperature and material phase may be found in [52, 42] . The function L is the laser function, defined in this work as a heat source, as done in [35, 34] .

The functions S p , S l , and S s are sigmoid functions selected to satisfy the following properties:

we may write the system in a slightly altered version of the form given by (1):

By summing (44)- (46) , one observes that:

at each (x, t). By initializing φ p = 1, φ l = 0, φ s = 0 (consistent with the physics), then expect that, at all (x, t):

It follows then that the system satisfies the hypotheses of Theorem 2; the phase fractions satisfy a mass conservation condition, while the temperature equation (43) , coupled to the phase-fraction dynamics, is separate from the compartmental system.

We consider a square domain assumed to be 400 microns in height and width, over a time interval of 7.5 milliseconds. We activate the laser function, defined as:

where (x c , y c ) are 200 and 400 microns, respectively. The laser is activated for the first 1.5 milliseconds of the time period, and then switched off for the remainder, for the cooling phase.

In space, we discretize with P 1 finite elements for the temperature and phase fields, over a mesh with size h = 8 microns. For the boundary conditions, we fix the temperature along the bottom of the domain as 293.15 K, and in time, we use an implicit BE scheme with ∆t = 0.0125 milliseconds. These spatial and temporal discretizations were selected after mesh convergence analyses showed that an adequate level of spatial and temporal convergence were reached. In Fig. 7 , we show the solutions of the system (43)- (46) at an instance in time. We see that the laser, defined at the top of the domain, provides a concentrated heat source, which melts the powder phase. At the instance shown, the laser has been deactivated, and we see the solidification process in action; while some areas of the melted powder remain hot enough such that the liquid phase dominates, the outer area has cooled enough to solidify. In Fig. 8 , we show this cooling/solidification process in more detail, showing the liquid (top) and solid (bottom) phases at four different instances in time, observing the change from liquid into solid.

We perform a similar analysis as in the previous example. We first compute the reference solution over the full time interval. Then, we define a training limit and reconstruct the DMD solutions for both the coupled and uncoupled DMD methods. In this instance, we defined the training limit to be t = 3.75 seconds, halfway through the simulation. As before, we are interested in seeing how the coupled and uncoupled DMD compare, in terms of both their error behavior and their ability to properly respect the physics (e.g., mass conservation).

In Fig. 9 we plot the error behavior in time for the temperature and phase fraction fields for both the coupled and uncoupled cases. For the temperature field, the error behavior for both reconstructions is very similar. This is unsurprising, as the temperature field is somewhat independent of the compartmental physics, and therefore can be reconstructed with reasonable accuracy without considering the phase fractions. However, it is worth noting that, when considered jointly with the phase fractions, the error behavior shows little difference. We note that, in both the coupled On the top row, we show the liquid fraction at four instances in time, with the bottom row depicting the solid fraction at the same four instances. As the material cools, the liquid begins to solidify, beginning first with the outer areas, before cooling and reaching a near-complete solid in the fourth instance shown. Figure 9 : Error behavior in time for the additive manufacturing simulation. We see very similar error performance for both the coupled and uncoupled DMD for the temperature field, suggesting a relative independence from the phase compartment fields. However, for the phase fractions, the uncoupled DMD completely fails, with only the powder compartment showing reasonable behavior. The liquid and solid compartments, perhaps due to being initialized at zero, fail for both reconstruction and prediction. However, when applying the coupled DMD reconstruction, this effect is not an issue.

and uncoupled case, the temperature error performance is very strong; it remains under 1% for a long period of time after the end of the training period, and is still under 10% at the end of the time interval.

When considering the error in the phase fraction compartments, we see a clear advantage when using the coupled approach. The error remains low for both the reconstruction and prediction period for each phase fraction for coupled DMD. We observe errors under 1% for the training period, and, as with the temperature field, they remain low throughout the time interval, even long after the training period is finished. However, for the uncoupled case, we observe failure in both the reconstruction and prediction. The solid and liquid phases, perhaps due to their initial conditions, remain near zero for the entirety of the time interval. In this sense, the long-term error behavior for the liquid compartment, while superficially appearing acceptable, is misleading; as the physical compartment tends to zero, the error appears acceptable. However, during the period in which the liquid compartment displays nonzero behavior, the uncoupled DMD fails to reproduce this dynamics, and the resulting error during this time is large. Figure 10 : Melt pool profile for the reference case and the two different DMD approaches. We see reasonable agreement with the reference for both DMD methods; however, the coupled DMD approach nonetheless gives a demonstrably better result, matching the reference nearly perfectly, while the uncoupled DMD solution deviates more from the reference.

We next examine the melt pool depth, often a quantity of much interest in simulations of additive manufacturing [52] . We measure this by observing the powder phase fraction along the vertical centerline for the reference, coupled, and uncoupled DMD solutions. We plot these results in Fig. 10 . In this case, both the coupled and uncoupled DMD provide reasonable solutions (as expected given the error behavior of the powder fraction in Fig. 9 ). However, the coupled approach nonetheless demonstrates superior performance, as the coupled approach shows noticeably more deviation from the reference solution towards the top of the domain.

In Fig. 11 , we finally examine the mass-conservation between the reference computation and the different DMD methods. As in the previous example, the coupled DMD reconstruction is mass-conservative throughout the considered interval, consistent with the theory. We note additionally that, as this model also incorporates a temperature field, independent of the compartment structure, this example is a numerical confirmation of Theorem 2. We again see that, when using an uncoupled DMD approach, the mass is not conserved.

In this work, we examined two approaches for the dynamic mode decomposition method on compartmental PDE systems: the uncoupled approach, in which DMD is applied to each compartment separately, and the coupled approach, in which DMD is applied to all compartments together. The need for this study was motivated by the inconsistent performance observed for DMD as a predictive tool in a number of situations. We established theoretically that if a compartmental system is mass-conservative, as is often the case, we may expect DMD performed on the coupled system to respect this property, even when extrapolated in time. It was also shown that if this coupling holds only for a subset of the considered variables, we still nonetheless expect this property to hold when DMD is applied to the entire system. Based on these results, the error behavior may be improved as long as the compartments of the reconstructed system remain non-negative. Figure 11 : Mass conservation for the additive manufacturing problem. The mass is properly conserved with the coupled DMD analysis, confirming our theory. However, as in the SIRD model, the uncoupled DMD does not conserve mass, and the mass loss is quite pronounced in both the reconstruction and prediction periods.

We proceeded to show, through numerical examples, that coupled DMD can provide accurate reconstructions, even in cases when uncoupled DMD fails to give useful results, thus confirming our theoretical results. For an epidemiological model of COVID-19, it was shown that DMD provides accurate short-term predictions up to several weeks after the end of the training window. Further, the theoretical result concerning mass conservation was confirmed, and the error behavior was significantly better when compared to the uncoupled case, as expected. We proceeded to look at another application, motivated by additive manufacturing, which considered the evolution of a temperature field and the resulting change of phase between the powder, liquid, and solid states. This further confirmed our theory, and the uncoupled DMD was able to show good reconstructive and predictive ability and maintain the important mass-conservation properties.

The dynamic mode decomposition and applied Koopman theory is a very active area of modern research, with important works being published at an incredible rate. Such works include both theoretical analysis of DMD, as well as practical new algorithms and applications. Improving on the work shown here can be extended in both directions. In the application and practical setting, such as high order DMD [38] , parameterized DMD, bagging for DMD [48] , multilevel DMD [17] , adaptive DMD [22] , and examining the performance of the coupled-vs.-uncoupled formulations is an important and interesting future path. Fast algorithms which may help alleviating possible additional costs of coupled DMD are also an important area of study; for systems with a very large number of compartments, the cost and memory requirements of coupled DMD may become prohibitive. Systematic ways to quantify the coupling such that one may reduce the cost, by, for instance, determining if one may obtain the advantages of coupled DMD by performing it only on a subset of compartments, are desired. On the theoretical side, there are several mathematical results that we believe can be proven rigorously. As shown here, if one can assume non-negativity of the compartments in the DMD reconstruction, then the error of these reconstructions may remain bounded. However, in general, the reconstruction may not be non-negative, and as such, a formulation of DMD that could guarantee this property would potentially prove a valuable tool in using DMD for longer-term time extrapolations.

A review of variational multiscale methods for the simulation of turbulent incompressible flows

Understanding mass transfer directions via data-driven models with application to mobile phone data

Dynamic mode decomposition for compressive system identification

Dynamic mode decomposition in adaptive mesh refinement and coarsening simulations

Computational fluid-structure interaction: methods and applications

Numerical methods for delay differential equations

Hyperbolic models for the spread of epidemics on networks: kinetic description and numerical methods

Least-squares finite element method for a meso-scale model of the spread of covid-19

Modeling and simulating the spatial spread of an epidemic through multiscale kinetic transport equations

Mathematical models in epidemiology

On the formulation of epidemic models (an appraisal of kermack and mckendrick)

Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations

Data-Driven Science and Engineering: Machine Learning, Dynamical Systems, and Control

Machine learning for fluid mechanics

Dynamic Mode Decomposition Analysis of High-Fidelity CFD Simulations of the Sinus Ventilation. Flow, Turbulence and Combustion

Variational multiscale methods in computational fluid dynamics

Dynamic mode decomposition for forecasting and analysis of power grid load data

The Approximation of One Matrix by Another of Lower Rank

Randomized matrix decompositions using R

Data-driven nonlinear aeroelastic models of morphing wings for control: Data-driven nonlinear aeroelastic models

Adaptive mesh refinement and coarsening for diffusion-reaction epidemiological models

Assessing the spatio-temporal spread of COVID-19 via compartmental models with diffusion in Italy, USA, and Brazil

Delay differential equations for the spatially-resolved simulation of epidemics with specific application to covid-19

A deep collocation method for the bending analysis of Kirchhoff plate

Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions

Scientific Computing: An Introductory Survey, Revised Second Edition

Multiscale and stabilized methods. Encyclopedia of Computational Mechanics Second Edition

Low-rank dynamic mode decomposition: Optimal solution in polynomial-time

Bayesian-based predictions of covid-19 evolution in texas using multispecies mixture-theoretic continuum models

Numerical simulation of a susceptible-exposed-infectious space-continuous model for the spread of rabies in raccoons across a realistic landscape

Containing papers of a mathematical and physical character

Discovering spatial-temporal patterns of covid-19 pandemic in south korea: A data-driven approach using dynamic mode decomposition

libmesh: a c++ library for parallel adaptive mesh refinement/coarsening simulations

Accurate prediction of melt pool shapes in laser powder bed fusion by the non-linear temperature equation including phase changes

A hierarchical computational model for moving thermal loads and phase changes with applications to selective laser melting

Dynamic mode decomposition: data-driven modeling of complex systems

Multiresolution dynamic mode decomposition

Higher order dynamic mode decomposition

Deep learning of material transport in complex neurite networks

Reaction diffusion system prediction based on convolutional neural network

Prediction accuracy of dynamic mode decomposition

Recommended Values of Thermophysical Properties for Selected Commercial Alloys

Mathematical Biology II: Spatial Models and Biomedical Application

Discovering dynamic patterns from infectious disease data using dynamic mode decomposition

Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations

Recent developments in variational multiscale methods for large-eddy simulation of turbulent flow

Spectral analysis of nonlinear flows

Bagging, optimized dynamic mode decomposition (bop-dmd) for robust, stable forecasting with spatial and temporal uncertainty-quantification

Dynamic mode decomposition of numerical and experimental data

Application of the dynamic mode decomposition to experimental data

Image-based modelling for Adolescent Idiopathic Scoliosis: Mechanistic machine learning analysis and prediction

Numerical solution of additive manufacturing problems using a two-level method

Simulating the spread of COVID-19 via a spatially-resolved susceptible-exposed-infected-recovered-deceased (SEIRD) model with heterogeneous diffusion

Diffusion-reaction compartmental models formulated in a continuum mechanics framework: application to covid-19, mathematical analysis, and numerical study

Dynamic mode decomposition (dmd) of solids

The authors declare that they have no conflict of interest.