Dynamic mode decomposition in adaptive mesh refinement and coarsening simulations.

Dynamic mode decomposition (DMD) is a powerful data-driven method used to extract spatio-temporal coherent structures that dictate a given dynamical system. The method consists of stacking collected temporal snapshots into a matrix and mapping the nonlinear dynamics using a linear operator. The classical procedure considers that snapshots possess the same dimensionality for all the observable data. However, this often does not occur in numerical simulations with adaptive mesh refinement/coarsening schemes (AMR/C). This paper proposes a strategy to enable DMD to extract features from observations with different mesh topologies and dimensions, such as those found in AMR/C simulations. For this purpose, the adaptive snapshots are projected onto the same reference function space, enabling the use of snapshot-based methods such as DMD. The present strategy is applied to challenging AMR/C simulations: a continuous diffusion-reaction epidemiological model for COVID-19, a density-driven gravity current simulation, and a bubble rising problem. We also evaluate the DMD efficiency to reconstruct the dynamics and some relevant quantities of interest. In particular, for the SEIRD model and the bubble rising problem, we evaluate DMD's ability to extrapolate in time (short-time future estimates).

Introduction

Data-driven methods are currently revolutionizing the modeling, prediction, and control of complex systems. Increasingly, researchers are considering data-driven approaches for a diverse range of complex systems, such as turbulent flows, climate sciences, epidemiology, finance, robotics, and many other different applications [12]. Even with the availability of better hardware and advances in techniques and algorithms, numerical simulations of these systems are still resource-demanding: strong nonlinearities, multiple scales, and large dimensionalities are typical examples of complexities found in modern applications. With the assembly of modern mathematical methods, unprecedented data availability, and increasing computational resources, previously complex, challenging problems can now be tackled within the new research field entitled scientific machine learning (SciML).

SciML is a core component of artificial intelligence and computational technology that can be trained, with scientific data, to augment or automate human skills [6]. This emerging research area aims at the opportunities and challenges in the context of complex applications across science and engineering, and other interdisciplinary fields. A wide range of SciML methods can be categorized regarding the type of information available, and their intended use [13]. For instance, many state-of-the-art methods and applications within the SciML scope can be found in [13, 36, 45, 46, 52, 60]. In this study, we focus on Dynamic Mode Decomposition (DMD), an unsupervised SciML method that can extract the most dynamically relevant low-rank structures from large-dimensional data observed in dynamical systems. DMD can be applied to both numerical [56] and experimental data [57].

The standard DMD procedure for numerical simulations consists of stacking the m snapshots (discrete solutions in space for a given time step) in columns to create a matrix and map the dynamics using a linear operator. The DMD procedure assumes that the m collected snapshots have spatial dimensionality n, where , such that the snapshots matrix has dimension . This occurs when snapshots are obtained from experimental data with sensors in fixed positions or numerical simulations considering a fixed mesh. However, in many situations, this is not always achievable. For numerical simulations using Adaptive Mesh Refinement/Coarsening (AMR/C), for instance, spatial adaptivity leads to solutions computed in meshes that constantly change in time. Adaptive meshes lead to a different number of nodes, nodal coordinates and numbering, and mesh topologies. Although there are strategies regarding missing spatial data or subsampling for the DMD [14, 41], specific treatments related to AMR/C approaches in the finite element context are not explored in the literature. In the present study, we develop a strategy to project all the snapshots of a given simulation with different dimensionality onto a reference target mesh with minor accuracy loss, enabling the use of any SVD-based data-driven technique such as DMD.

This paper is structured as follows: Sect. 2 describes the relation between the discretization of PDEs in space and time and a dynamical system. This section introduces Dynamic Mode Decomposition, our method of choice for short-time future estimates and extrapolation. Section 3 describes our strategy to deal with simulations that consider AMR/C in their evolution (e.g., in the case of dimensionality of the output vector, as well as of mesh topology and/or node numbering, varying in time). In Sect. 4 we describe the numerical applications in this study: the use of DMD on a continuous SEIRD model for COVID-19 and two fluid dynamics problems, a density-driven gravity current, and a bubble rising problem. We show efficiency and accuracy results for the signal reconstruction. Moreover, future time step predictions using DMD are evaluated for the SEIRD model and the bubble rising problem. In Sect. 5, we draw our final remarks and conclusions.

Numerical methods and dynamic mode decomposition

Solving partial differential equations (PDEs) using fast, accurate, reliable, and robust methods is crucial for many industrial and scientific applications. Several methods (such as finite elements, finite differences, and many others) are responsible for approximating the infinite-dimensional PDEs into finite-dimensional spaces. The discretization of these equations allows the process to be automated. In the present study, we focus on using the finite element method for spatial discretization of the PDEs. That is, consider a generic transient parametric PDE such aswhere is a nonlinear operator, is a vector of parameters (e.g., diffusion, density, viscosity, etc.), is a given function and the solution is a function of spatial coordinates , temporal coordinates t, and parameters such that . The equation is equipped with boundary and initial conditionswhere and are the Dirichlet and Neumann boundary conditions, respectively, T is the final time, and is the initial condition for . The domain is bounded by Lipschitz continuous boundaries , and is the unit outward normal to . The union of boundaries and domain is represented as . The standard finite element method consists of discretizing into a mesh composed of nodes and elements. Each element has its domain and boundary . The weak form of the system can be obtained by integrating Eq. (1) in its strong form against a weighting function , where is the Sobolev space of the square-integrable functions with an integrable first weak derivative, and applying the divergence theorem. Being the space of polynomials of degree equal or less than k over , the function spaces are defined asTherefore, the semi-discrete finite element formulation for Eq. (1) is: find such that :where the inner product over the domain is indicated by . The weak form given by Eq. (5) naturally accommodates several finite element formulations, from Galerkin to Variational Multiscale methods [1, 8, 19, 39, 53]. After the temporal discretization of Eq. (5), the equation can be translated into a discrete-time dynamic system. In this system, the state vector at the time instant can be written such that:where represents the discrete-time flow map of the system and incorporates information regarding the parameters , mesh size, solver tolerances, etc. In the present study, we consider that the measurements of the system are the state vectors themselves, that is, . Analyzing the evolution in time of a discretized PDE as a dynamical system is a key concept for introducing Dynamic Mode Decomposition (DMD).

DMD is an equation-free, data-driven method that provides accurate assessments of the dominant structures in a given complex system [43]. DMD provides a decomposition of data into spatio-temporal modes that correlates the data across spatial features and also associates them to unique temporal Fourier modes. The main idea of the method is to efficiently compute the regression of linear/nonlinear terms to a least-square linear dynamics approximation from experimental or numerical observable data. Despite its first appearance in the fluid dynamics context [55, 56], DMD has been used in many other applications such as epidemiology [51], biomechanics [16], urban mobility [3], climate [44] and aeroelasticity [28], especially in structure extraction from data and control-oriented methods.

We can now apply DMD on the dynamical system described in Eq. (6). Consider a dataset containing the observations in time of the dynamical system for , where is the total number of observations. The datasetcan be split into two datasets and . DMD consists on finding the best fit approximation of the linear mapping that transforms dataset into dataset , that is,The computation of can be done as , where is the Moore–Penrose pseudoinverse of . However, we avoid the computation of the full matrix since is a 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 aswhere and are the left and right singular vectors and is a diagonal matrix with real, non-negative, and decreasing entries named singular values. The singular values are hierarchical and can be interpreted in terms of how much the singular vectors influence the original matrix . For the DMD procedure, considering the Eckart–Young Theorem [23], the optimal low-rank update approximation matrix , when subjected to a truncation rank r, can be written aswhere is a matrix containing the first r columns of , contains the first r columns of , and is the diagonal matrix containing the first r singular values. The pseudoinverse can be approximated asand, instead of computing , we can obtain , a projection of asNote that is unitarily similar to . 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 [40]. Now we can compute the eigendecomposition of :where is a diagonal matrix containing the discrete eigenvalues and the matrix contains the eigenvectors of . The DMD basis can be written asand the signal reconstruction asbeing the vector containing the projected initial conditions such that , and is a diagonal matrix whose entries are the continuous eigenvalues , where is the time step size between the outputs. In the present study, , where , being the time step size used in the temporal integration of the PDEs. For instance, if one chooses to output the solutions once every two time steps, the time step size between the two observations will be two times larger than the time step size used to compute time integration, that is, . It is important to mention that the snapshots sampling frequency affects directly the DMD’s ability to capture the dynamics. For lower dominant frequencies, a larger is more adequate, while smaller values of are required for capturing rapid dynamics [56]. The form of (15) can be regarded as a generalization of the Sturm–Liouville expansion for a differential problem:where and are the ith Sturm–Liouville eigenfunctions and eigenvalues for a given differential operator.

DMD can be seen as a dimensionality reduction method due to its inherent ability to extract the r most relevant dynamical modes, where r is often much smaller than the snapshot matrix rank m. However, a strategy to determine the number of relevant modes is not straightforward, and is an active topic in DMD research [59]. Even though the choice of r for DMD may require some trial and error, some techniques can be used to find a good starting point. A hard threshold technique [43] consists of choosing r such thatwhere is a tolerance threshold, set, e.g., to . This method implies that more than of the variance in the data is retained by the approximation. For the case where DMD is used on experimental (or numerical but noisy) data, more sophisticated solutions are presented in the literature [22, 29].

Another important consideration when using DMD is the SVD algorithm. The SVD can represent a significant part of the computational effort, meaning that improvements in the SVD performance lead to significant CPU time gains. For many SVD-based methods (such as DMD), there is no need to compute the SVD for the whole matrix, since the method aims at extracting the first r dominant structures in the matrix. Many algorithms are designed in this direction to make this computation more efficient. One important contribution in this direction is seen in [58], where the method of snapshots was proposed, paving the way to more algorithms and techniques. In this paper, we employ the randomized SVD (rSVD) algorithm [25, 38], a non-deterministic algorithm able to compute the near-optimal low-rank approximation of a given large dataset with good efficiency.

DMD on adapted meshes

Regarding numerical methods to approximate PDEs, the use of finer meshes in finite element discretizations usually leads to more accurate solutions. On the other hand, reducing the number of equations in the nonlinear system is crucial for efficiency, especially considering that the optimal computational complexity of a single physics transient finite element simulation is [15]. The duality between the two statements describes a well-known trade-off between accuracy and efficiency in the finite element context. Despite a considerable research effort in the past decades, strategies to generate tailored meshes to maximize the accuracy while minimizing the computational effort are still an open research topic. Milestones addressing this subject are finite element a posteriori error estimators/indicators [2], techniques such as adaptation, interpolation [17, 47] and projection [18, 27, 50], as well as libraries and frameworks containing automated versions of AMR/C techniques [4, 42].

The general structure of the AMR/C scheme is given in Algorithm 1 and illustrated in Fig. 1. Three criteria are fundamental in an AMR/C algorithm: remeshing, flagging, and stopping. The remeshing criterion defines whether the computed solution at a given time step requires remeshing driven by global a posteriori error estimators (or indicators) and/or by calling the refinement/coarsening procedure at every j time steps. Next, all mesh elements are visited and flagged for refinement or coarsening. The element flagging criterion is often represented by local a posteriori error estimators or indicators, i.e. flux-jumps of the solution gradient. In possession of the flagged elements, the remeshing algorithm is invoked. Finally, the previous mesh solution must be projected or interpolated into the target mesh. This is done by the projection/interpolation algorithm. The whole process is repeated until the stopping criterion is achieved. This criterion could be error equidistribution (until a certain threshold), a given element size, the maximum number of elements in the mesh, or the number of refinement levels. Note that the meshes generated in Algorithm 1, guided by the error estimation procedure, have different dimensions (number of degrees of freedom) and different topologies, which poses difficulties for DMD (or, in fact, for any snapshot-based method).

Fig. 1

Illustration of a mesh refinement procedure. A local a posteriori error estimator or indicator flags an element for refinement (in green) using the solution computed in the mesh on the left. The mesh is refined (or coarsened) according to the flagged elements, and the process can be restarted until a given criterion is met (error level, element size, maximum number of elements, etc.). Note that the initial and the final mesh differ in the number of degrees of freedom and topology

The structure of the exact DMD algorithm relies on the fact that the measurements of the state vectors, i.e., have the same dimensionality. This structure could be exemplified in numerical experiments as fixed discretizations in space, i.e., fixed meshes or static sensors in experimental data. However, finite element simulations equipped with AMR/C strategies provide solutions in different function spaces, depending on the mesh used to compute the solution on a given time step. Spatial adaptivity on transient finite element simulations leads to meshes with a different number of nodes, numbering, and nodal coordinates. It can also lead to different mesh topologies and structures. For that reasons, snapshots obtained by AMR/C simulations cannot be stacked in columns to construct the snapshot matrix. Even if one considers an AMR/C strategy that restricts the adaptive meshes to preserve the dimensionality of the snapshots, the difference between the nodal coordinates of the various meshes will lead to misleading dynamics captured by DMD. In the present study, we circumvent this issue by considering the -projection [50] of the numerical simulation results for different meshes into a reference target mesh. Recent work on reducing these projection costs for related reduced-order modeling techniques (though not DMD) may be found in [26, 32, 37]. The projection or interpolation of numerical solutions between finite element meshes is a well-known computational mechanics subject. Many issues regarding boundary conditions, data visualization, or coupling arise from this kind of problem. The choice of a proper method to successfully project functions in different finite element spaces is not a trivial task since conservation may not be satisfied [18, 27, 50]. This issue is not addressed in the present study since the mesh projection occurs as a post-processing phase after computing the AMR/C solutions. Therefore, it cannot lead to cumulative errors. Furthermore, we are also careful to choose reference target meshes with characteristic length equivalent with the existent in the AMR/C meshes to avoid any major accuracy losses. For the sake of generality, we consider the -projection as our method of choice, and it can be defined as follows [61]. Assuming that the solution on the donor mesh to be projected onto the target mesh must satisfy the orthogonality condition,where is a finite-dimensional subspace of defined by the target mesh and the interpolant is the optimal interpolant in the -norm for . The orthogonal projection can be defined in terms of the following linear systemwhere is the mass matrix and is the projection matrix. The mass matrix is usual in finite element computations. The matrix, however, can present technical difficulties (see “Appendix”).

Therefore, our strategy consists of applying the -projection onto a tailored reference target mesh capable of representing the many scales in time and space of all the snapshots. This routine is inserted in the code and invoked after the adaptive procedure in Algorithm 1 at every output time step . This tailored reference mesh is described in this work as target reference mesh and should not be confused with the target mesh generated during the AMR/C procedure. The new solution (computed on the adaptive donor mesh) is projected onto the reference target mesh and exported as a simulation output. In this work, we export as output files for visualization purposes, although there is no restriction on stacking the snapshots on the snapshot matrix during the simulation runtime or dumping only on disk. We considered Gmsh [30], an open-source robust mesh generator, as our software of choice for defining and creating the target meshes for this study. The output for each time step is the snapshot with constant dimensionality n such that all nodes in space are correctly mapped and capturing the dynamics existent in the system. The procedure is summarized in Algorithm 2. This strategy is relatively simple since the -projection consists of solving a linear system where the generated matrix is a mass matrix and requires no extra outputs for storing the projected solutions since the projection can be applied right after the AMR/C code. The mass matrix is generated in the finite element context by a self-adjoint operator, enabling more efficient solvers. In terms of versatility, the -projection method is flexible because a solution obtained for a given mesh can be naturally projected onto reference target meshes with different topologies and dimensionalities. Also, since the mesh projection is a vital part of AMR/C algorithms, finite element algorithms frequently present efficient implementations of interpolation or projection techniques. Figure 2 shows an example where a solution obtained by an adaptive mesh simulation is projected onto two meshes with different topologies.

Fig. 2

Comparison of different -projection examples on structured and unstructured meshes

For this example, we consider the square domain and the function u defined in such thatThis function is approximated on a structured finite element mesh discretized into cells where each cell is divided into two triangular elements. An AMR/C procedure is invoked to refine three times the transition between and , creating a new mesh containing 1672 elements and 857 nodes. Figure 2a shows the new mesh generated after the AMR/C procedure and the function approximated by the resulting function space. Two meshes—one structured and another unstructured—are considered for projection. The structured mesh contains cells, resulting in 12,800 triangular elements and 6561 nodes. The element sizes of the structured mesh are similar to the smallest elements in the adaptive mesh. The unstructured mesh presents a smaller characteristic length than the structured mesh and contains 17,088 elements and 8705 nodes. Figure 2b, c shows the proposed meshes and the projections of the solution onto the new finite element spaces. Note that the projected solutions are fairly accurate since their infinity norm is in good agreement with the adaptive mesh solution’s infinity norm.

Remark

In this study, the projection is carried out inside the finite element simulation since the projection computational cost is practically negligible in comparison with the overall time required for solving nonlinear systems of a complex numerical simulation. However, if one does not have access to the finite element simulation codes used, the reconstruction of the solution of each time step can still be done off-line. Output files of various formats contain information regarding the mesh used (such as nodal coordinates and connectivities) for visualization purposes. By properly reading these files, the solutions can be reconstructed under a finite element framework (such as FEniCS [4] or libMesh [42]) or on a code developed by the user. However, this totally non-intrusive approach can significantly increase the computational cost due to several I/O operations that are often extremely low compared to computational intensive operations. Since the mesh constantly varies in an AMR/C finite element simulation, each solution obtained in the simulation has to be imported, reconstructed under its original mesh, and projected. After the projection, the user can choose to dump the solution in the disk or stack the snapshots on the snapshot matrix. For the first case, another I/O operation would be invoked. This non-intrusive strategy is especially suitable (and restricted) to be used with DMD, which is also a non-intrusive algorithm. For intrusive snapshot-based methods such as Proper Orthogonal Decomposition (POD), one needs to project finite element matrices onto the computed basis and, therefore, requires access to the code [62].

Numerical experiments

In this section, for the sake of generality, we apply our method in several applications, with different systems of equations, numerical formulations, mesh topologies, spatial dimensions, refinement criteria, and finite element libraries. We compare the results between the AMR/C solution, the fixed mesh solution, and the DMD results for all cases. First, we test the DMD short-time future prediction capabilities on a continuous SEIRD model for COVID-19, a nonlinear system of diffusion–reaction equations. The equations are considered using a Galerkin finite element discretization and are solved using libMesh [42], a high-performance C++ finite element library. The error estimators, refinement/coarsening strategies built-in on libMesh can be seen on [34, 63]. Also, the -projection algorithm is embedded in libMesh. We explore the results in one and two spatial dimensions, where the 1D case is a hypothetical example, and the 2D case describes the COVID-19 evolution in the Lombardy region in Italy [34, 63, 64]. The simulation obtains the results for 44 and 60 days for the 1D and 2D cases, respectively. Since we want to predict 14 days in the future, we only feed DMD with snapshot matrices containing the first 30 and 46 days in the 1D and 2D cases, respectively, such that the predicted results can be compared with the results obtained in the simulations. To apply DMD to the SEIRD data obtained on an AMR/C simulation, we project the adaptive solution onto a mesh with characteristic length compatible with the obtained simulation results. The discretization of the reference mesh used for the projection is as fine as the final refinement level on the adaptive meshes. Next, we consider two fluid dynamics applications where the governing equations for both cases are advection-dominated. To circumvent the LBB condition and spurious oscillations regarding dominant advection, these equations consider the residual-based variational multiscale (RBVMS) formulation [1, 8, 19, 35, 39, 53] on a finite element discretization. We consider the use of DMD on a 2D density-driven gravity flow and a 3D bubble rising simulation. The density-driven gravity current is modeled by the coupling of the incompressible Navier–Stokes equation and the advection–diffusion equation. We consider a lock-exchange problem, where a tank is filled with two fluids of different densities, separated by a lock. The simulation starts when the lock is removed, and the difference in the densities of the fluids generates the driving forces responsible for the motion of the fluids. Unlike the other examples in this work, the implementation of this numerical test is made on the FEniCS v.2019.1 framework [4], a high-performance Python/C++ finite element library. The refinement/coarsening algorithm and -projection algorithm for this example are part of the framework. We consider an interface-tracking error indicator for the AMR/C simulations, and the mesh is refined following the bisection method [54]. DMD is considered to reconstruct the solution, and the results are compared to the fixed mesh results and the results obtained by AMR/C simulations. Finally, we extend our analysis to a 3D bubble rising case [33], a two-phase incompressible flow problem where the interface is captured by the convected level-set method. This model is implemented on libMesh, taking advantage of the same refinement/coarsening strategies as well as the projection algorithm used in the SEIRD numerical tests. In this example, we test the projection of the adaptive solutions onto three different meshes containing the different scales existent on the AMR/C simulation and evaluate the results. We consider DMD to predict the bubble geometry and dynamics for a short time in the future.

For all the numerical tests proposed, we evaluate the results in terms of efficiency and accuracy. For efficiency purposes, we compute the ratio between the computational time required to run the finite element simulations and the time required to run DMD separately. We refer to this quotient as speedup. The finite element code is responsible for computing the snapshots and projecting the results onto the reference target meshes proposed, while the DMD code imports the output files, extract the snapshots, computes the approximation, and outputs the results. Also, we provide a table describing how the projection routine affects the overall computational time required for the simulations for all examples. In terms of accuracy, we evaluate the results in terms of overall relative error where is the snapshots matrix, is a matrix comprising the approximations obtained by DMD and denotes the Frobenius norm. A more detailed analysis is also done in terms of relative error in time . For that, we plot the curves of the relative errors (in terms of -norm) of each snapshot. We compute for snapshots, where denotes the -norm. Also, to avoid unphysical results, some quantities of interest are evaluated and compared using both simulation and approximation results. For the SEIRD model, we plot the results regarding the total population. This quantity should be constant in time according to the hypothesis of the model. For the lock-exchange simulation, we compute the mass during the simulation and the front position. Since the simulation is considered on a closed tank, the mass must be kept constant during the simulation. For the 3D bubble rising problem, we plot the quantities of interest related to the geometry (volume and sphericity) and dynamics (center of mass and rise velocity) of the bubble. In this case, we test meshes with different minimum characteristic lengths. The results are shown and discussed below.

Continuous SEIRD model for COVID-19

The outbreak of COVID-19 in 2020 has led to a surge in interest in the mathematical modeling of infectious diseases. This new virus is responsible for infecting millions of people worldwide and impacting the economy in an unprecedented way. Therefore, numerical simulation of the virus’ dynamics may help provide short-term prediction models for forecasting the number of future cases. In this perspective, it is possible to develop strategic planning in the public health system to avoid deaths and manage patients.

Disease transmission may be modeled as 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 [11]. Here, we work with a spatio-temporal SEIRD model, presented in [34, 63, 64], and given by,where , , , , and denote the densities of the susceptible, exposed, infected, recovered, and deceased populations, respectively. The sum of all the compartments with the exception of is represented by which is the total living population. characterizes the Allee effect (persons), that takes into account the tendency of outbreaks to cluster around large populations, and denote the transmission rates between symptomatic and susceptible individuals and asymptomatic and susceptible individuals, respectively (units days), denotes the incubation period (units days), corresponds to the asymptomatic recovery rate (units days), the symptomatic recovery rate (units days), represents the mortality rate (units days), and , , , are the diffusion parameters of the different population groups as denoted by the sub-scripted letters (units km persons days). Note that all these parameters can be considered time and space-dependent. We also compute the compartment c, the cumulative field of the i compartment.

For the numerical solution of (21)–(25), 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 Differentiation Formula (BDF2), which offers second-order accuracy while remaining unconditionally stable. We implement the whole model in libMesh [42]. We additionally make use of AMR/C, allowing us to resolve multiple scales. One may find more details about the methods in [34, 63].

Reproducing a 1D model

First, we use a simple 1D continuous SEIRD model for COVID-19 with adaptive mesh refinement to validate the -projection and DMD. This example was introduced in [63] and reproduced in [34]. Basically, we consider a 1D region with initial conditions that represents a large population centered around with no exposed persons and a small population centered around with some exposed individuals, as shown in Fig. 3. Thus, we set and as follows:We further set , , and . We also enforce homogeneous Neumann boundary conditions at and a zero-population Dirichlet boundary condition at for all model compartments. The latter represents a non-populated area at .

Fig. 3

Initial conditions for the 1D model

Following [34, 63] we set days, dayspersons, days, days and days, , , , and kmpersonsdays. The time step size is defined as days and we consider the mesh projection and outputs at every time step, that is, .

We use an adapted mesh with initially 125 elements, and after the refinement, the smallest element has a size 0.002. At the beginning of the simulation, we refine uniformly the whole domain into two levels and, after that, we apply the adaptive mesh refinement every 4 time steps. The idea is that the AMR/C strategy will keep this spatial resolution on more dynamically relevant regions while coarsening other regions in the domain. As a target reference mesh, we consider the uniformly refined mesh, such that all the domain contains the minimal spatial resolution obtained by the AMR/C simulation.

The results of this simulation are validated against the results from [34, 63]. Figure 4 show the solution of the 1D SEIRD example at days for the fixed mesh simulation, the adaptive mesh simulation, and the projection of the adaptive solution onto the target reference mesh used in the fixed mesh simulation. We observe a good agreement between the solutions. This agreement is a positive indicator that the DMD can be used on the projected meshes with little to no error compared to the DMD on a fixed mesh simulation. In terms of efficiency, Table 1 shows the computational effort required for the reference mesh projection embedded on the adaptive finite element code in comparison with the time used for the simulation code itself (AMR/C FEM).

Fig. 4

Solution at days for the fixed mesh solution, AMR solution and the respective projection onto a reference mesh for the 1D SEIRD example. The reference mesh was built with characteristic length similar to the smaller elements in the adaptive mesh

Table 1

Absolute and relative time required for the projection routine in comparison with the adaptive finite element simulation code (AMR/C FEM) for the SEIRD model in the 1D case

Code part	Absolute time (s)	Relative time (%)
AMR/C FEM	1262.89	98.53
Mesh projection	18.88	1.47

After some initial experiments, we set for all compartments. In this example, the finite element simulation computes the results for 44 days, that is, 176 time steps. However, we only consider the first 30 days in the snapshot matrix for reconstruction. DMD is set to approximate the results for a further 14 days so the results can be compared with the 44 days results obtained in the finite element simulation. In other words, we want to evaluate the DMD ability to predict the COVID-19 scenario two weeks in the future given the data of the last 30 days. For numerical reasons, we considered an initial 3 days shift in the snapshots since some compartments are initialized with zeros, affecting how DMD captures the dynamics. The 44th day solution for the adaptive simulation and the 44th day prediction considering DMD are seen in Fig. 5. We observe good agreement between predictions and the numerical solutions of the simulations for most compartments. The speedup and overall relative error between the DMD approximation and the snapshots are seen in Table 2.

Fig. 5

Solution at days for the AMR simulation solution and the 14 days projection using DMD for the 1D SEIRD example

Table 2

Relative error between reconstructed (and predicted) data and the projected snapshots

Compartments	Relative error	Speedup
s	1.590×10^-3	1.079×10^3
e	2.574×10^-2	1.316×10^3
i	1.162×10^-2	1.225×10^3
r	1.439×10^-2	1.264×10^3
d	2.001×10^-2	1.331×10^3
c	1.286×10^-2	1.341×10^3

The Lombardy region

We extend our analysis by solving the continuous SEIRD model and applying DMD to a 2D real world domain that is the Lombardy region in Italy. The spread of the COVID-19 has been studied in this region using the continuous SEIRD model with accurate results [63, 64]. Here, we reproduce this simulation with the solver developed in [34] which invokes adaptive mesh refinement every 4 time steps. We use the same parameters as shown in [63, 64]. It is important to point that, in this simulation, the transmission rates and diffusion parameters vary with time in order to reproduce the effects of restrictions during the simulated period.

For this simulation, an unstructured mesh is considered due to the complex geometry imposed by the domain. The mesh is generated using Gmsh and is uniformly refined as the simulation starts. After refining the whole mesh in one level, the mesh presents a minimum spatial resolution of approximately 1 km. This procedure allows the solver to coarsen the regions where no significant dynamics are observed while preserving the scales of the regions of interest. Figure 6 shows the variation of the number of nodes in time for the AMR/C strategy. The coarsening approach improves the simulation performance significantly since the average number of nodes (and respectively, the number of equations ) used in the adaptive simulation is approximately half with respect to the case of a fixed mesh with the same spatial resolution considered for the entire domain. In this example, the reference target mesh is the uniformly refined unstructured mesh considered in the early stages of the simulation, presenting 13,158 nodes and 25,340 elements. The simulation considers a time step size of days for the numerical integration and days for the observations. Initial conditions for the Lombardy domain are the same presented on [63, 64] and are seen in Fig. 7, while compartments r and d are initialized to zero.

Fig. 6

Number of mesh nodes in time for the adaptive solution and the proposed reference mesh for the 2D SEIRD example

Fig. 7

Initial conditions for the SEIRD model in the Lombardy case

We then proceed to run both adaptive and fixed mesh simulations and, to apply DMD in the adaptive mesh results, we consider the proposed projection scheme. Figure 8 shows the s compartment solution at days for both simulations and the projected adaptive solution onto the reference mesh, revealing that the results are in good agreement. In terms of efficiency, Table 3 shows results for the computational time required for the projection compared to the finite element code. That said, the adaptive snapshots can now be assembled into a snapshot matrix for the DMD reconstruction and prediction.

Fig. 8

Solution for the susceptible compartment at days obtained using an adaptive mesh and its respective projection onto a fixed reference mesh

Table 3

Absolute and relative time required for the projection routine in comparison with the adaptive finite element simulation code (AMR/C FEM) for the SEIRD model in the Lombardy case

Code part	Absolute time (s)	Relative time (%)
AMR/C FEM	2107.24	98.52
Mesh projection	31.60	1.48

The DMD analysis is made in the same way as presented in the 1D case: the simulation outputs the projected snapshots for the first 60 days (240 snapshots). The snapshot matrix assembles the information regarding 46 days of simulations, while DMD approximates the results for 60 days. The idea is to predict two weeks in the future, given the data observed in the past 46 days. For this example, we consider the SVD truncation at for all compartments. Again, an initial 3 days shift in the data is considered to avoid issues with the compartments initialized to zero. Results for the 60th day comparing the computed numerical solutions and the DMD predictions are seen in Figs. 9 and 10. We present the projected solutions, the DMD prediction, and the relative error in space between the two results from left to right. We can note that most compartments show results in agreement with the simulations, while the exposed compartment reveals more pronounced differences than the other compartments.

Fig. 9

Comparison between computed and predicted solutions at days for the susceptible, exposed, and infected compartments

Fig. 10

Comparison between computed and predicted solutions at days for the recovered, deceased, and cumulative infected compartments

Figure 11 shows the relative error in time between the DMD results and the projected snapshots. The first thing we notice is that the curves are different for each compartment. This discrepancy occurs due to the different parameters for each equation in the SEIRD model, which largely affects the dynamics of the system. The dynamics for each compartment are different since each compartment presents different coupling, diffusion, and reaction parameters. Also, regarding this issue, since the parameters are time and space-dependent, sudden changes in their values can affect the dynamics of the system as well as DMD’s dynamics mapping ability. Some sudden changes in the s and e compartments related to stricter public policies considered to reduce the transmission rates (parameters and ) are incorporated into the model. Since the variation in the parameters is not introduced smoothly, DMD’s ability to map sudden changes in the dynamics of the system is reflected by the existence of some spikes on the curves of the relative errors in time. Comparing the reconstruction and prediction stages, we observe that the errors tend to grow as soon as the prediction stage starts (dashed line). The exposed compartment, which yielded most of the oscillations due to parameter changing on the reconstruction stage, presented the same behavior on the prediction phase around day 49. We also note that the exposed compartment yields a larger relative error for the 60th day in comparison with the other compartments. Table 4 shows the overall relative error and the speedups for the six compartments approximations. Comparing these results with the results presented in Figs. 9, 10 and 11, we can conclude that the predictions are reasonably accurate in comparison with the numerical solutions, specially when considering the time required for calculation.

Fig. 11

Relative error for all compartments between numerical simulation snapshots and DMD reconstruction and prediction. The dashed line represents the beginning of the DMD prediction stage

Table 4

Relative error between reconstructed (and predicted) data and the computed snapshots and speedup between DMD and the numerical simulation

Compartments	Relative error (η_F)	Speedup
s	1.345×10^-3	822.61
e	5.187×10^-2	938.09
i	1.304×10^-2	755.97
r	7.316×10^-3	1036.93
d	1.097×10^-2	977.91
c	1.251×10^-2	977.03

Another important analysis to be done is the conservation property of the continuous SEIRD model. As mentioned before, the standard -projection does not guarantee conservation among the projections. Figure 12 shows the total population during the simulation, normalized by the total population modeled in the initial conditions. The total population is computed as the sum of the integral of the compartments (excluding c) divided by the sum of the integral of the elements of the mesh. Since the SEIRD model does not consider any population growth, the value must be theoretically constant for all the simulations. From the figure, we observe that the population is kept constant during all the adaptive simulation, and this was preserved by the projected solutions and the DMD reconstruction stage. That is, we can note that the -projection does not yield conservation issues in this example. For the prediction phase, DMD preserves the total population for several days in the future. However, it presents a slight increase (around ) for predictions over 10 days, which does not affect the results significantly. This increase can be explained by the relative errors behavior, observed in Fig. 11, as DMD computes future estimates.

Fig. 12

Population conservation for both adaptive and projected results

Fluid dynamics

This section evaluates the DMD use on two cases involving AMR/C in computational fluid dynamics: the reconstruction of a 2D density-driven gravity current simulation and the temporal prediction on a 3D rising bubble. Different from the previous cases, the test cases presented in this section are advection-dominant. To approximate the governing equations, we use a finite element RBVMS formulation [1, 19, 35, 39, 53]. In the first case, we reconstruct the solution and evaluate important quantities of interest regarding density-driven gravity flows. On the bubble rising simulation, we show how our strategy works on a 3D mesh, and we evaluate the DMD ability to predict the quantities of interest of the rising bubble in time.

Density-driven gravity flow

In this section, we consider a long numerical simulation that consists of a lock-exchange between two fluids, the heavy fluid, A, and the lighter fluid, B, based on the numerical example in [49]. The difference between their densities is such that the Boussinesq hypothesis is considered valid. Moreover, particles in the heavy fluid have negligible inertia and are much smaller than the smallest length scales of the buoyancy-induced fluid motion. Thus, the dimensionless governing equations arewhere is the fluid velocity, is the concentration field, p is the pressure, is the vector pointing in the direction of gravity, is the Schmidt number and is the Grashof number, two dimensionless numbers that relate viscous effects with diffusion and buoyancy effects, respectively. A Grashof number of this magnitude indicates a turbulent flow. The field is the concentration and is responsible for mapping the evolution of fluid interactions. The time step size considered for this simulation consists on s for a total simulation time of s with an output frequency of s. We consider a tank, that is, a rectangular domain with length m, height m. The boundary conditions for this case are no-slip for the velocity and no-flux for the transport equation, and the initial conditions are such that the heavy fluid is represented as a column with dimensions m 2m located at the left border of the tank and the light fluid fills the rest of the domain. Figure 13 illustrates the domain and the initial conditions.

Fig. 13

Scheme illustrating the initial conditions for the density-driven gravity flow example

To solve the governing equations, we implement the RBVMS formulation [35] for Eq. (28) using the FEniCS 2019.1 [4] framework to generate the snapshots for this example. The adaptive mesh refinement procedure returns the solution for and . For this example, we only consider the snapshots of for our calculations, that is, an overall data reduction of . Details of the formulation of the problem can be found in [7, 35]. We consider a fixed mesh simulation with nodes and cells, where each cell is divided into two linear triangles. We consider an interface-tracking adaptive mesh error indicator that flags and refines the mesh where the two fluids interact for the adaptive mesh simulation. The error indicator for the mesh refinement is being larger than a given tolerance. For this purpose, a mesh containing cells is considered, the interfaces are refined considering two levels of refinement, and the mesh refinement is invoked at every time step. Figure 14 shows the results at s for both fixed and adaptive mesh simulations and the projection of the adaptive solution onto the same mesh used in the fixed mesh simulation.

Fig. 14

Results and mesh for the first 8 m of the domain at s

AMR/C in this problem is advisable since the dynamics are predominant on the interface between the fluids. Most of the domain is not affected in the early stages of the simulation, and the use of fine meshes outside these regions may represent unnecessary computational effort. Figure 15 shows the number of nodes in the mesh during the simulation time for the adaptive mesh compared with the fixed nodes mesh.

Fig. 15

Number of mesh nodes in time for the adaptive solution and the proposed reference mesh for the density-driven gravity current example

We observe that the number of nodes in the adaptive simulation is smaller than the fixed mesh for all the simulation time. Table 5 shows the time required for the simulation to run in comparison with the time spent on projecting the solutions onto the reference target mesh.

Table 5

Absolute and relative time required for the projection routine in comparison with the adaptive finite element simulation code (AMR/C FEM) for the lock-exchange example

Code part	Absolute time (s)	Relative time (%)
AMR/C FEM	15748.11	94.94
Mesh projection	839.42	5.06

Now we proceed applying DMD to the projected solution and reconstructing the solutions. Figure 16 shows the relative error for the reconstruction using different values of the rank r. The results are also confirmed in Table 6, where the relative error between the reconstructed and projected snapshot matrices is shown and the speedup computed for each case. We observe that, for increasing values of r, the relative error decreases for all steps and affects the overall relative error of the matrices. That is, inserting more structures in the DMD basis yields better accuracy in terms of overall relative error. As for the speedup, we notice that it decreases for increasing values of r. Such a decrease occurs because larger values of r directly affect the rSVD algorithm performance [7, 25] and increase the dimensions of the matrices for the computation of the DMD basis.

Fig. 16

Relative error for the reconstruction considering different values of the rank r

Table 6

Relative error between reconstructed data and the snapshots and speedup between DMD and the numerical simulation

Rank r	Relative error (η_F)	Speedup
50	1.081×10^-1	558.53
100	3.472×10^-2	445.91
150	6.207×10^-3	230.67
200	2.024×10^-3	204.74
250	9.604×10^-4	189.17

Most importantly, we can evaluate the quantities of interest common to density-driven gravity currents. Figure 17 shows the front position and mass conservation for both simulations and the respective reconstructions. We note that, for these quantities all reconstructions are in extremely good agreement with those computed with the fixed and adaptive meshes, approaching them, as expected, for higher values r.

Fig. 17

Front position and mass conservation for the fixed mesh and adaptive mesh simulations and reconstructions with the target mesh

Bubble rising problem

We now study a bubble rising 3D benchmark, whose task is to track the evolution of a three-dimensional bubble rising in a liquid column. The initial configuration is described in Fig. 18.

Fig. 18

Initial configuration and boundary conditions for the bubble rising problem

For this problem, we couple the Navier–Stokes equations with an interface capturing method called convected level-set [20, 33, 65]. The convected level-set is a method used to represent the interface between two phases and, by a convection equation, to move the interface as the flow evolves. A force that has an important role in bubble problems is the surface tension , which is applied using the Continuum Surface Model (CSF) [10].

We write the governing equations in their dimensional form aswhere is the density, is the dynamic viscosity, is the acceleration of gravity vector, is the level-set function, is a penalty constant, , and S a function related to the level-set signed distance function.

For the temporal integration, we apply the Backward–Euler method to the Navier–Stokes equations, while for the convected level-set, we use the BDF2 method. One may find more details about the governing equations and methods in [33].

The initial configuration consists of a spherical bubble of radius m centered at [0.5, 0.5, 0.5] m in a m domain. The no-slip boundary condition is applied to all boundaries. Table 7 lists the parameters used for this simulation.

Table 7

Rising bubble data

Computational domain	1×1×2	(m)
Grid sizes	0.100 to 0.025	(m)
Number of time steps	240	(–)
Time step	0.0125	s
Bubble radius	0.25	m
Initial bubble position	(x,y)=(0.5,0.5,0.5)	m
Liquid density	1000	kg/m^3
Liquid viscosity	10	kg/(ms)
Gas density	100	kg/m^3
Gas viscosity	1	kg/(ms)
Surface tension	24.5	N/m
Gravity	0.98	m/s^2

We use an adapted mesh, initially with cells, with each cell divided into 6 linear tetrahedra. We refine the initial region where the bubble is located into two levels and, after the refinement, the smallest element has a size of 0.025 m. The adaptive mesh refinement is based on the flux jump of the level-set function error, in which . We apply the adaptive mesh refinement every four time steps. The interface is modeled with , and the time step size is defined as s. We output the projected solutions at every 2 time steps such that s. In this example, we consider three tetrahedral meshes for our projection strategy, presented in Fig. 19. The three meshes named coarse, intermediate and fine, present characteristic lengths similar to the three scales existent in the refinement levels of the adaptive mesh simulation. The coarse mesh represents the initial mesh on the adaptive simulation, with 12,000 elements and 2541 nodes. The intermediate mesh presents smaller elements equivalent to the generated elements after the first refinement level on the AMR/C simulation, totalizing 96,000 elements and 18,081 nodes and the fine mesh contains the smallest scales presented on the adaptive numerical solution with 768,000 elements and 136,161 nodes. The figures show the level-set solution on half of the domain and at s for the adaptive mesh solution and the respective projections. We also present on Table 8 the time required for the projection of the solutions of on the three meshes in terms of absolute and relative time. We can see that the projection time is very small.

Fig. 19

Level-set solution detail at s and projection to the coarse, medium and fine meshes

Table 8

Absolute and relative time required for the projection routine in comparison with the adaptive finite element simulation code (AMR/C FEM) for the bubble rising example

Target mesh	Code part	Absolute time (s)	Relative time (%)
Coarse	AMR/C FEM	17345.92	99.97
	Mesh projection	4.41	0.03
Intermediate	AMR/C FEM	17390.03	99.85
	Mesh projection	24.56	0.14
Fine	AMR/C FEM	17653.70	98.91
	Mesh projection	194.87	1.09

To verify if the bubble geometry and dynamics are being preserved during the mesh projections, we evaluate the quantities of interest such as the bubble volume, sphericity, center of mass, and rise velocity. Bubble volume and sphericity are related to the bubble geometry, while the position of the bubble center of mass and rise velocity regard the bubble dynamics. The results are compared for the AMR/C simulation output and the projections on Fig. 20. The quantities of interest are computed using the adaptive snapshots as well as the projected snapshots and compared with each other. Regarding the geometry quantities of interest, we observe that the coarse mesh projection solution affects the bubble’s geometry, leading to bad results regarding the bubble volume and sphericity compared with the AMR/C simulation. For the intermediate and the fine meshes, the bubble’s geometry is not largely affected, and the results are compatible with the AMR/C results. As for the quantities of interest related to the bubble dynamics, no significant difference is observed on the projection of the three meshes. We observe that the center of mass is not affected by the projections compared with the simulation results. As for the rise velocity, we observe some minor differences in all projection cases.

Fig. 20

Comparison between the simulation and projection of the 3D rising bubble quantities of interest

The simulation proceeds until s, yielding a dataset containing 240 snapshots regarding the solution of for each target mesh from the projections. We do not consider the use of velocities and pressure in the DMD analysis, reducing the required data in . We consider the results for the first 2.75s to construct the basis and predict the last 0.25 s. We compare the DMD results in terms of relative error between the snapshot matrix and the obtained solutions for each case. The DMD solution for the coarse mesh is compared to the projected adaptive solution onto the coarse mesh and so forth. Results are evaluated for multiple values of r, such that . The results regarding accuracy and performance for multiples values of r and the three meshes are presented in Table 9.

Table 9

Relative error between reconstructed data and the projected snapshots and speedup between DMD and the numerical simulation

Rank r	Mesh	Rel. error (η_F)	Speedup
5	Coarse	6.222×10^-2	1.33×10^5
	Intermediate	6.655×10^-2	7.95×10^4
	Fine	6.865×10^-2	1.53×10^4
10	Coarse	2.384×10^-2	1.21×10^5
	Intermediate	2.461×10^-2	3.60×10^4
	Fine	2.490×10^-2	1.22×10^4
15	Coarse	1.415×10^-2	1.15×10^5
	Intermediate	1.402×10^-2	2.46×10^4
	Fine	1.429×10^-2	7.55×10^3
30	Coarse	8.900×10^-3	1.09×10^5
	Intermediate	8.125×10^-3	3.29×10^4
	Fine	7.687×10^-3	6.66×10^3
45	Coarse	6.382×10^-3	1.02×10^5
	Intermediate	5.820×10^-3	5.35×10^4
	Fine	5.776×10^-3	5.80×10^3
60	Coarse	5.659×10^-3	8.74×10^4
	Intermediate	5.444×10^-3	1.45×10^4
	Fine	5.518×10^-3	5.53×10^3

Results presented for multiple values of r

The results for are shown in Fig. 21. We observe that the errors are stable for the reconstruction case, that is, the DMD solution before s. From that point, shown as a dashed line on the figure, the errors begin to grow exponentially for each predicted time step, while still remaining below 1% until around 2.9 s. We observe that the errors in the reconstruction case are different regarding the mesh used. That is, the errors are larger with respect to the minimum characteristic length of the projection meshes. However, we observe that the errors grow at the same rate on the prediction phase independently of spatial discretization.

Fig. 21

Relative error for the rising bubble example for the coarse, intermediate, and fine mesh solutions. The dashed line defines the start of the prediction phase

We also show the results considering for reconstruction at s and prediction at s in Fig. 22. We compare the DMD results for the three projected meshes, the adaptive mesh solution, and a fixed mesh solution. The fixed mesh solution is obtained by running the simulation with a mesh of 768,000 elements and 136,161 nodes, as the fine mesh used in the projection. We observe initially that the bubble geometry is better defined on the reconstruction than on the prediction figure. This better definition is directly related to the errors observed in Fig. 21. We observe that the coarse mesh results do not capture the bubble geometry with the same accuracy as the intermediate and fine meshes for the reconstruction results. As for the intermediate and fine meshes, they present similar results in comparison with the projected solutions. However, when we observe the prediction figure, we observe that instabilities inherent to DMD arise on the bubble contour, affecting the bubble geometry for the intermediate and fine mesh.

Fig. 22

Bubble contour at the vertical mid plane for the signal reconstruction ( s) and prediction ( s) last steps

We now proceed comparing the results in terms of quantities of interest for the DMD results. Figure 23 shows the bubble volume, sphericity, center of mass, and rise velocity for the DMD results compared to the adaptive solution results. The same issue regarding sphericity on the coarse mesh projection is observed on the coarse mesh DMD results. However, for the intermediate and fine meshes, the values match the results observed for the projection. We observe results in conformity for the center of mass and rise velocity as well. In terms of prediction, we observe that DMD accurately predicts the volume and the center of mass evolution. As for the other quantities of interest, the increasing exponential errors in the DMD prediction structure affect the quantities of interest for long time future predictions.

Fig. 23

Comparison between the simulation and DMD signal plus prediction of the 3D rising bubble quantities of interest. The dashed line marks the beginning of the prediction regime for the DMD

Conclusions

In this work, we propose a strategy to enable data-driven snapshot-based methods on finite element solutions obtained by AMR/C simulations. The use of AMR/C algorithms in finite element approximations of PDEs is known to reduce memory usage and to increase the efficiency of the simulations without compromising the accuracy. The simulation adapts the mesh with the evolution of the solution, refining regions of interest and coarsening regions that are not of interest. The adaptation process leads to different meshes during the simulations, and the solutions have different dimensions and topologies (or different connectivities and nodal coordinates) whenever the AMR/C algorithm is invoked. Snapshots with different dimensions prevent the use of snapshot-based algorithms such as DMD, where the snapshot matrix is built by stacking the solutions for each observation in columns. In this study, we considered a strategy to project the adaptive solutions on reference meshes such that all snapshots present the same dimensions and nodal indices, enabling the construction of snapshot matrices. The method employed to project the solutions onto the target reference mesh is the -projection, a simple, fast and versatile approach. The -projection is a common strategy present in several finite element libraries and frameworks and consists basically of solving a linear system, where the matrix is obtained by a self-adjoint operator enabling the use of efficient solvers. Despite presenting drawbacks, especially regarding properties’ conservation, we investigate the use of the -projection as a postprocessing tool without any relevant issues. By postprocessing, we mean that the -projection algorithm is invoked only to insert the solutions on a reference function space for each time step and output the files. This strategy does not yield significant additional computational effort as we observe that the -projection routine required around of the computational time required for the adaptive finite element code to run in most cases and in the worst case. When the source code is not available to invoke the -projection, one can construct the solutions from the output files and project them onto the target reference mesh in a complete non-intrusive workflow.

We test the algorithm on several models presenting different dynamics and underlying physics. First, we present the results for DMD on a continuous SEIRD model for COVID-19 for fictional data in 1D and real data for Lombardy, Italy, in 2D. The idea of considering short-time future estimates on COVID-19 models could improve the decision-making of public policies to avoid further contamination, and the use of AMR/C in the simulations, coupled with the presented projection scheme and DMD, can lead to fast and reliable predictions. The simulations are implemented using the libMesh library using its refinement built-in functions. We compare the solution, its projection in the reference mesh, and the adaptive solution and notice that no relevant errors are found for the projection strategy. In this example, the -projection consists of and of the total time required for the code to run for the 1D and 2D cases, respectively. As for the use of short-time future estimates, we consider DMD for predicting two weeks in the future given a set of snapshots. For the 1D case, we feed DMD with snapshots covering 30 days, while for the 2D case, the snapshot matrix comprises 46 days of simulation results. The DMD results are presented in terms of efficiency and accuracy compared to the projected solutions, and we observe that the predictions are mostly in good agreement except for the exposed compartment. To test the approach’s agnosticism regarding the dynamics of the systems, we consider two different fluid dynamics problems: a 2D density-driven gravity current and a 3D bubble rising problem. We use a 3000-snapshot simulation to generate a basis containing the spatio-temporal coherent structures for the density-driven gravity flow. The lock-exchange simulation is solved using the FEniCS framework. The computational effort regarding the projection, in this case, corresponds to of the overall computational time. As for the DMD analysis, we evaluate the results for different values of the rank r and observe that, as we increase it, in this case, the overall relative errors decrease significantly. However, when evaluating the quantities of interest such as mass conservation and front position, even the worst case yields good approximations. For the 3D example, we use the libMesh library, testing three different reference meshes considering the three element sizes existent in the adaptive simulation. We investigate the occurrence of projecting the adaptive solution onto meshes that do not necessarily guarantee that all the scales obtained in the solution are preserved in terms of quantities of interest related to the bubble geometry and dynamics. We observe that the use of -projection on coarser meshes leads to issues regarding the shape of the bubble, affecting quantities of interest such as bubble volume and sphericity. When considering the -projection on meshes with approximately the same spatial resolution as the finest scale in the adaptive meshes, these effects are mitigated. As for the bubble dynamics, even coarser meshes reveal similar results for the bubble center of mass and rise velocity in comparison with finer mesh projections and the adaptive solution itself. In terms of computational effort, the -projection corresponds to , and of the total simulation time for the coarse, intermediate, and fine meshes, respectively. Proceeding to the DMD analysis, we present the results for a 220 time step reconstructions and a 20 time step prediction for the bubble rising problem. We test the approximations for multiple values of r being our best approximations. The results computed by DMD are presented in terms of efficiency and accuracy compared to the projected solutions. We also reevaluate the quantities of interest, now with the DMD approximation. For the reconstruction, results are practically identical to those observed for the projection, while for the prediction, the center of mass and bubble volume yield good results while sphericity and rise velocity are not in agreement with those obtained in the simulation.

Summarizing our findings, we highlight that the present approach, despite its simplicity:The results observed in this study reveal that the use of mesh projections for adaptive solutions preserves the dynamics inherent to the finite element solutions. However, this approach on large problems could be inefficient or prohibitive since the mesh projection on a sufficiently fine target reference mesh can lead to unfeasible DMD storage and computation times. As a future work, aiming to circumvent this issue, we plan to explore the use of data compression on DMD, once used for streaming purposes [24], in the mesh projection context. Also, it is known that DMD approximates the nonlinear dynamics existent in the dataset using a linear operator, being fragile to problems presenting a strong nonlinearity [43]. To deal with this, several Koopman-based frameworks, such as the Extended DMD [66], the Kernel DMD [67], and LANDO [5], use nonlinear observables to improve temporal predictions. In the future, we plan to evaluate our projection approach on nonlinear observables for further investigations.

is versatile since every finite element library that presents AMR/C is often equipped with efficient projection algorithms;

is fast in the sense that the -projection only requires solving a sparse, well-conditioned linear system where is the mass matrix, obtained by a self-adjoint operator, enabling the use of efficient solvers;

is framework-agnostic, as noted in the results of our simulations implemented on libMesh and FEniCS, requiring only the standard procedures already in place in both libraries;

is problem-agnostic, as observed on meshes with different spatial dimensions (1D, 2D, and 3D), topologies (structured and unstructured), and equations (convection-dominated and reaction–diffusion systems);

preserves the dynamics computed on the adaptive snapshots for all the tested examples.