# Parametric Nonlinear Model Reduction UsingK-Means Clustering for Miscible Flow Simulation

Authors: Norapon Sukuntee; Saifon Chaturantabut
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

This work considers the model order reduction approach for parametrized viscous fingering in a horizontal flow through a 2D porous media domain. A technique for constructing an optimal low-dimensional basis for a multidimensional parameter domain is introduced by combiningK-means clustering with proper orthogonal decomposition (POD). In particular, we first randomly generate parameter vectors in multidimensional parameter domain of interest. Next, we perform the K-means clustering algorithm on these parameter vectors to find the centroids. POD basis is then generated from the solutions of the parametrized systems corresponding to these parameter centroids. The resulting POD basis is then used with Galerkin projection to construct reduced-order systems for various parameter vectors in the given domain together with applying the discrete empirical interpolation method (DEIM) to further reduce the computational complexity in nonlinear terms of the miscible flow model. The numerical results with varying different parameters are demonstrated to be efficient in decreasing simulation time while maintaining accuracy compared to the full-order model for various parameter values.


---

## Body

## 1. Introduction

Viscous fingering is the formation of patterns in a morphologically unstable interface between two fluids in a porous medium. It generally refers to the onset and evolution of instabilities that occur in the displacement of fluids [1]. The mathematical description of viscous fingering can be given in terms of several partial differential equations. Many numerical methods are used to find the approximation of pressure, velocity, concentration, or temperature. In [2], the coupled set of partial differential equations is approximated by the semidiscrete SUPG stabilized finite element formulation plus a discontinuity capturing technique to improve stability around the moving sharp fronts. The pressure equation is discretized by the standard Galerkin method, and a postprocessing scheme is used to improve the numerical evaluation of Darcy’s velocity. Branston et al. developed a finite volume numerical simulator with a level of certainty where fingering perturbations in stochastic fields do not result into severe channeling as discussed in [3]. The simulator development involves discretization of convection–diffusion and pressure equations with total variation diminishing and central differencing schemes, respectively. Correlated stochastic fields for simulation runs were created using convolution with a Gaussian filter achieved through 2D fast and inverse Fourier transform, Dykstra–Parsons coefficient, random porosity generator, and autocorrelation lengths. In [4], the Oldroyd model has been studied using a generalized incompressible Navier-Stokes equation via the interpolating stabilized element free Galerkin technique. In this work, the Crank-Nicolson method is first utilized to discrete the temporal direction, and then, the meshless element free Galerkin technique is employed to approximate the spatial direction. The equation contains the integral term which is approximated by using a difference scheme. Although the numerical methods discussed earlier are effective and accurate, the resulting simulations generally require enormous memory storage and consume excessively computational time. In order to manage this problem, model order reduction techniques have been proposed in many previous works for various applications.Proper orthogonal decomposition (POD) is a commonly used technique for model order reduction. POD is used to generate the optimal basis that can describe the dominant features of a dynamical system. This basis is utilized with the Galerkin projection to construct a reduced-order model. As in [5], a surrogate model based on the proper orthogonal decomposition is developed in order to enable fast and reliable evaluations of aerodynamic fields. In 2017, Chang et al. introduced a reduced-order model for wall shear stress in abdominal aortic aneurysms by proper orthogonal decomposition [6]. A more interesting detail on POD can be found in [7–9].Using POD still has a limitation when the model has complicated nonlinear property. POD cannot handle the complexity issue arising in computing nonlinear terms. A powerful technique called the discrete empirical interpolation method (DEIM) is used here to overcome this problem. In particular, the DEIM is employed to construct selected interpolation indices, which will allow much less expensive evaluation for the nonlinear terms. The DEIM can be used to combine with POD to improve the efficiency of the nonlinear reduced-order model. In [10], the POD-DEIM approach has been applied to the groundwater equation and the time-dependent incompressible Navier-Strokes equation with variable density based on the meshfree RBF-FD technique. In [11], POD-DEIM is applied to a dynamic 2D catalytic reactor model. There are also many other applications of POD-DEIM on various problems, which can be found, for example, in [12–14].In this paper, we apply the POD-DEIM technique on a finite difference discretization of a nonlinear miscible flow with viscous fingering in porous media. This model contains many parameters which can be changed in many different situations. Our additional effort is to construct a reduced-order model, which can be used to approximate the solution with any given parameter in the domain of interest. This issue in model reduction is crucial for various parametrized dynamical systems. There are a number of works that have been done for parametric model order reduction (PMOR). In [15], the PMOR method for laminated composite structures is presented. The parametrization is carried out for the angle between the layers of structure and the thickness of the structure. In [16], this work proposed the PMOR scheme for second-order systems that does not require a priori sampling of the parameter space. For the flow model, PMOR for probabilistic analysis of unsteady aerodynamic applications is discussed in [17]. Also in [18], the PMOR approach is applied to the fluid flow governed by the Navier-Stokes equations at varying Reynolds number. In this work, we use K-means clustering to efficiently construct a reduced system for a parametrized miscible flow model.This paper is arranged as follows. We first formulate the model equations in Section2. In Section 3, a general setup for constructing POD-DEIM reduced systems is discussed. Section 4 introduces a parametrized nonlinear model reduction using K-means with the POD-DEIM approach for the miscible flow model. In Section 5, the numerical experiments are used to illustrate the efficiency of the proposed approach. We finally conclude the results and discuss some future research directions in Section 6.


## 2. Problem Formulation

We consider a viscous fingering in a horizontal flow through a 2D porous media domain [19] with a constant permeability K given by
(1)∇·u=1,(2)∇P=−μKu,(3)∂c∂t+u·∇c=D∇2c+fc,(4)ρcp∂T∂t+u·∇T=DT∇2T+−ΔHfc,where u=u1,u2∈ℝ2 is the velocity with components in x and y coordinates, P is the global pressure, c is the concentration of the injected fluid, T is the temperature, μ is the viscosity depending on c and T, fc is the reaction rate, and ρ, cp, D, DT, and ΔH are assumed to be constant. To solve (1)–(4), a general procedure nondimensionalizes the system which will be converted into the equations in terms of stream function ψ and vorticity ω as follows:
(5)∇2ψ=−ω,(6)ω=−Rcψxcx+ψycy+cy+RTψxTx+ψyTy+Ty,(7)∂c∂t+ψycx−ψxcy=∇2c+Dafc,(8)∂T∂t+ψyTx−ψxTy=Le∇2T+sgnϕDafc,where Rc and RT are constant parameters determining the effects of concentration and temperature to the viscosity, Da (Damköhler number) and Le (Lewis number) are dimensionless parameters,
(9)sgnϕ=1;exothermic,−1;endothermic,and fc=−cc−1c+d. These equations (5)–(8) are defined on the spatial domain x,y∈Ω=0,αPe×0,Pe⊂ℝ2, where Pe (Péclet number) is a dimensionless constant aspect ratio and time domain t∈0,tf⊂ℝ. For this system, the boundary conditions are
(10)ψx,0,t=ψx,Pe,t,ψyx,0,t=ψyx,Pe,t,cx,0,t=cx,Pe,t,cyx,0,t=cyx,Pe,t,Tx,0,t=Tx,Pe,t,Tyx,0,t=Tyx,Pe,t,along top and bottom sides. The boundary conditions are
(11)ψ0,y,t=0,ψαPe,y,t=0,c0,y,t=1,cαPe,y,t=0,T0,y,t=1;exothermic,0;endothermic,TαPe,y,t=0;exothermic,1;endothermic,along the remaining two sides. The initial conditions are
(12)cx,y,0=Tx,y,0=1;x≤x^0;x>x^forally∈0,Pe,where x^=αPe/2 is the interface location. Let 0=x0<x1<⋯<xnx<xnx+1=αPe and 0=y0<y1<⋯<yny<yny+1=Pe be equally spaced grid points on x-axis and y-axis with Δx=αPe/nx+1 and Δy=Pe/ny+1, respectively. Define ψt, ωt, ct, and Tt as a vector of dimension n=nxny containing approximate solutions ψxi,yj,t, ωxi,yj,t, cxi,yj,t, and Txi,yj,t, respectively, for i=1,2,⋯,nx and j=1,2,⋯,ny. The central finite difference is applied to construct the spatial discretization of (5)–(8), and it becomes a system of nonlinear ODEs in the form of
(13)Mpdztdt=Apzt+dp+Ep;zt,t,where
(14)Mp=OOOOOOOOOOInOOOOIn,Ap=−A−InOOO−In−RcAyRTAyOOAOOOOLeA∈ℝ4n×4nare parametrized constant matrices;
(15)dp=OObLeb∈ℝ4nis a parametrized constant vector; Ep;·,·:ℝ4n×ℝ⟶ℝ4n is a nonlinear function with
(16)Ep;zt,t=O−RcNψt,ct+RTNψt,Tt−Fψt,ct+Dafct−Fψt,Tt+signϕDafct;zt=ψt;ωt;ct;Tt∈ℝ4nis the vector of state variables, and p=Rc;RT;Da;Le;Pe is the vector of parameters. Note that A,Ax,Ay∈ℝn×n are constant coefficient matrices, O∈ℝn×n is the zero matrix, In∈ℝn×n is the identity matrix, b, bx∈ℝn are vectors indicating the boundary conditions, and N,F are nonlinear functions that maps some state variables to ℝn. We strongly recommend that the reader should read more details in the derivation of (13) in Appendix A. The solution z will generally depend on the parameter vector p, and therefore, we will also denote zt as zp;t=ψp;t;ωp;t;cp;t;Tp;t. More details on discretization and numerical schemes can be found in [19, 20].


## 3. Reduced-Order Modeling by the POD-DEIM Approach

This section provides a brief overview of the general procedure for constructing a reduced-order model of a differential equation using proper orthogonal decomposition (POD) and a discrete empirical interpolation method (DEIM).


### 3.1. Subspace from Proper Orthogonal Decomposition

Consider a matrix of snapshotsZ=z1,z2,⋯,zns∈ℝn×ns where zj∈ℝn, j=1,2,⋯,ns. Suppose we want to approximate a snapshot zj by using a set of orthonormal vectors ϕ1,ϕ2,⋯,ϕk⊂ℝn, which has rank k<n. Then, the approximation is in the form
(17)zj≈Φkq,where Φk=ϕ1,ϕ2,⋯,ϕk∈ℝn×k is a matrix of basis vectors and q=q1,q2,⋯,qkT∈ℝk is the vector of unknown coefficients. To find q, we use the fact that Φk has orthonormal columns, i.e., ΦkTΦk=Ik, and the minimum error occurs when the residual is orthogonal to the column span of Φk, i.e., ΦkTzj−Φkq=0, which implies that q=ΦTzj and the approximation becomes
(18)zj≈ΦkΦkTzj.Proper orthogonal decomposition (POD) provides an orthonormal basis that minimizes this approximation error in 2-norm for a given basis rankk≤k∗≔rankZ; i.e., POD basis is the optimal solution to the minimization problem
(19)minΦk∈ℝn×k∑j=1nszj−ΦkΦkTzj22suchthatΦkTΦk=Ik,where Ik∈ℝk×k is the identity matrix. It can be shown [21] that the POD basis defined above can be obtained from the left singular vector of the snapshot matrix Z. Let Z=UΣVT be the singular value decomposition of Z, where matrices U=u1,⋯,uk∗∈ℝn×k∗ and V=v1,⋯,vk∗∈ℝns×k∗ are matrices with orthonormal columns and Σ=diagσ1,⋯,σk∗∈ℝk∗×k∗ is a diagonal matrix with σ1≥σ2≥⋯≥σk∗>0. Then, the POD basis matrix of dimension k is Uk=u1,⋯,uk∈ℝn×k, k≤k∗. Moreover, this minimum error is given by [21] the sum of neglected singular values squared, i.e.,
(20)minΦk∈ℝn×k∑j=1nszj−ΦkΦkTzj22=∑j=1nszj−UkUkTzj22=∑ℓ=k+1k∗σℓ2.As a result, the decay of singular values can be used to specify the dimensionk of the subspace of the POD basis used in the approximation. In most existing applications, POD can provide a basis that explains the overall dynamics of a given snapshot set and we can use reduced dimension k that is much less than the original dimension n, while maintaining high accuracy in the approximation.


### 3.2. Discrete Empirical Interpolation Method (DEIM)

In the context of nonlinear model reduction for dynamical systems, the DEIM [22] is generally used for reducing the complexity for evaluating the projected nonlinear term in the POD-Galerkin reduced system. To illustrate this issue, we consider the projected nonlinear function E~·,·:ℝk×ℝ⟶ℝk of E·,·:ℝn×ℝ⟶ℝn defined by
(21)E~z~t,t≔UkTEUkz~t,t,z~:ℝ⟶ℝk.For convenience, the notation will be simplified by leaving the dependence on time parametert. Notice that, for a given z~, the computational complexity for computing the above nonlinear term depends on the dimension n when multiplying UkT with E. This can be avoided by using the DEIM approximation of the form
(22)E~z~≈UkTU¯rPTU¯r−1PTE,where the matrix U¯r is the POD basis matrix of rank r of the nonlinear snapshot matrix Emat=E1,E2,⋯,Ens, for Ej=Ezj,tj with zj coming from the available snapshot in ℝn (zj≈Ukz~ from POD approximation), j=1,⋯ns, and the matrix P is constructed from the interpolation indices ℘1,℘2,⋯,℘r, which are generated by the DEIM algorithm shown in Algorithm 1. In particular, P=e℘1,e℘1,⋯,e℘r∈ℝn×r, where e℘i=0,⋯,0,1,0,⋯,0T is the ℘i column of the identity matrix In∈ℝn×n for all i=1,2,⋯,r for r≤RankEmat. Note that the complexity for evaluating (22) can depend only on the reduced-order dimension r and k because the term UkTU¯rPTU¯r−1 can be precomputed for any given input vector z and the term PTE can be considered as selecting r components corresponding to the indices ℘1,℘2,⋯,℘r.Algorithm 1: DEIM.
Require: Inputs (i) dimension r<RankEmat, (ii) basis set u¯j=1r=u¯1,u¯2,⋯,u¯r⊂ℝn.Ensure: Outputs P∈ℝn×r, ℘→r=℘1,℘2,⋯,℘rT∈ℝr1:∣ρ∣,℘1=max∣u¯1∣2:U¯=u¯1, P¯=e℘1, ℘→=℘13:forj=2:rdo4:c=P¯TU−1P¯Tu¯j5:r¯=u¯j−u¯c6:∣ρ∣,℘j=max∣r∣7:U¯⟵U¯u¯j, P¯⟵P¯e℘j, ℘→⟵℘→℘j8:end for9:P=P¯:,1:r, ℘→r=℘→1:rFrom Algorithm1, the procedure initializes the first interpolation index ℘1∈1,2,⋯,n by using the first input basis u¯1 entry which has the largest magnitude. The rest of the indices ℘j for j=2,3,⋯,r are selected from the residual r¯=v¯j−Uc with the largest magnitude. More details on this procedure as well as the proof for invertibility of PTUr can be found in [22].


### 3.3. The POD-DEIM Reduced System

Consider an initial value problem of the form(23)Mdzdt=Az+d+Ez,zt0=z0,where M∈ℝn×n, A∈ℝn×n, and b∈ℝn are constant matrices; Ez is the nonlinear function; and z=zt∈ℝn is the state variable. The POD reduced system can be obtained by applying the Galerkin projection with the POD basis Uk of dimension k (see Section 3.1), and this system can be written in the form
(24)M~dz~dt=A~z~+d~+E~z~,z~t0=UkTz0,where M~=UkTMUk∈ℝk×k, A~=UkTAUk∈ℝk×k, and d~=UkTd∈ℝk are constant matrices; E~z~=UkTEUkz~ is the projected nonlinear function; and z~=z~t∈ℝk is the state variable that can be used to approximate the original state variable z by z≈Ukz~. Finally, after approximating the nonlinear term as shown in Section 3.2, we obtain the POD-DEIM reduced system in the following form:
(25)M~dz~dt=A~z~+d~+QE~pz~,z~t0=UkTz0,where Q=UkTU¯rPTU¯r−1 and E~pz~=PTEUkz~ as discussed in (22). The next section extends the previous setup for reduced-order modeling for the parametrized miscible flow model.


## 3.1. Subspace from Proper Orthogonal Decomposition

Consider a matrix of snapshotsZ=z1,z2,⋯,zns∈ℝn×ns where zj∈ℝn, j=1,2,⋯,ns. Suppose we want to approximate a snapshot zj by using a set of orthonormal vectors ϕ1,ϕ2,⋯,ϕk⊂ℝn, which has rank k<n. Then, the approximation is in the form
(17)zj≈Φkq,where Φk=ϕ1,ϕ2,⋯,ϕk∈ℝn×k is a matrix of basis vectors and q=q1,q2,⋯,qkT∈ℝk is the vector of unknown coefficients. To find q, we use the fact that Φk has orthonormal columns, i.e., ΦkTΦk=Ik, and the minimum error occurs when the residual is orthogonal to the column span of Φk, i.e., ΦkTzj−Φkq=0, which implies that q=ΦTzj and the approximation becomes
(18)zj≈ΦkΦkTzj.Proper orthogonal decomposition (POD) provides an orthonormal basis that minimizes this approximation error in 2-norm for a given basis rankk≤k∗≔rankZ; i.e., POD basis is the optimal solution to the minimization problem
(19)minΦk∈ℝn×k∑j=1nszj−ΦkΦkTzj22suchthatΦkTΦk=Ik,where Ik∈ℝk×k is the identity matrix. It can be shown [21] that the POD basis defined above can be obtained from the left singular vector of the snapshot matrix Z. Let Z=UΣVT be the singular value decomposition of Z, where matrices U=u1,⋯,uk∗∈ℝn×k∗ and V=v1,⋯,vk∗∈ℝns×k∗ are matrices with orthonormal columns and Σ=diagσ1,⋯,σk∗∈ℝk∗×k∗ is a diagonal matrix with σ1≥σ2≥⋯≥σk∗>0. Then, the POD basis matrix of dimension k is Uk=u1,⋯,uk∈ℝn×k, k≤k∗. Moreover, this minimum error is given by [21] the sum of neglected singular values squared, i.e.,
(20)minΦk∈ℝn×k∑j=1nszj−ΦkΦkTzj22=∑j=1nszj−UkUkTzj22=∑ℓ=k+1k∗σℓ2.As a result, the decay of singular values can be used to specify the dimensionk of the subspace of the POD basis used in the approximation. In most existing applications, POD can provide a basis that explains the overall dynamics of a given snapshot set and we can use reduced dimension k that is much less than the original dimension n, while maintaining high accuracy in the approximation.


## 3.2. Discrete Empirical Interpolation Method (DEIM)

In the context of nonlinear model reduction for dynamical systems, the DEIM [22] is generally used for reducing the complexity for evaluating the projected nonlinear term in the POD-Galerkin reduced system. To illustrate this issue, we consider the projected nonlinear function E~·,·:ℝk×ℝ⟶ℝk of E·,·:ℝn×ℝ⟶ℝn defined by
(21)E~z~t,t≔UkTEUkz~t,t,z~:ℝ⟶ℝk.For convenience, the notation will be simplified by leaving the dependence on time parametert. Notice that, for a given z~, the computational complexity for computing the above nonlinear term depends on the dimension n when multiplying UkT with E. This can be avoided by using the DEIM approximation of the form
(22)E~z~≈UkTU¯rPTU¯r−1PTE,where the matrix U¯r is the POD basis matrix of rank r of the nonlinear snapshot matrix Emat=E1,E2,⋯,Ens, for Ej=Ezj,tj with zj coming from the available snapshot in ℝn (zj≈Ukz~ from POD approximation), j=1,⋯ns, and the matrix P is constructed from the interpolation indices ℘1,℘2,⋯,℘r, which are generated by the DEIM algorithm shown in Algorithm 1. In particular, P=e℘1,e℘1,⋯,e℘r∈ℝn×r, where e℘i=0,⋯,0,1,0,⋯,0T is the ℘i column of the identity matrix In∈ℝn×n for all i=1,2,⋯,r for r≤RankEmat. Note that the complexity for evaluating (22) can depend only on the reduced-order dimension r and k because the term UkTU¯rPTU¯r−1 can be precomputed for any given input vector z and the term PTE can be considered as selecting r components corresponding to the indices ℘1,℘2,⋯,℘r.Algorithm 1: DEIM.
Require: Inputs (i) dimension r<RankEmat, (ii) basis set u¯j=1r=u¯1,u¯2,⋯,u¯r⊂ℝn.Ensure: Outputs P∈ℝn×r, ℘→r=℘1,℘2,⋯,℘rT∈ℝr1:∣ρ∣,℘1=max∣u¯1∣2:U¯=u¯1, P¯=e℘1, ℘→=℘13:forj=2:rdo4:c=P¯TU−1P¯Tu¯j5:r¯=u¯j−u¯c6:∣ρ∣,℘j=max∣r∣7:U¯⟵U¯u¯j, P¯⟵P¯e℘j, ℘→⟵℘→℘j8:end for9:P=P¯:,1:r, ℘→r=℘→1:rFrom Algorithm1, the procedure initializes the first interpolation index ℘1∈1,2,⋯,n by using the first input basis u¯1 entry which has the largest magnitude. The rest of the indices ℘j for j=2,3,⋯,r are selected from the residual r¯=v¯j−Uc with the largest magnitude. More details on this procedure as well as the proof for invertibility of PTUr can be found in [22].


## 3.3. The POD-DEIM Reduced System

Consider an initial value problem of the form(23)Mdzdt=Az+d+Ez,zt0=z0,where M∈ℝn×n, A∈ℝn×n, and b∈ℝn are constant matrices; Ez is the nonlinear function; and z=zt∈ℝn is the state variable. The POD reduced system can be obtained by applying the Galerkin projection with the POD basis Uk of dimension k (see Section 3.1), and this system can be written in the form
(24)M~dz~dt=A~z~+d~+E~z~,z~t0=UkTz0,where M~=UkTMUk∈ℝk×k, A~=UkTAUk∈ℝk×k, and d~=UkTd∈ℝk are constant matrices; E~z~=UkTEUkz~ is the projected nonlinear function; and z~=z~t∈ℝk is the state variable that can be used to approximate the original state variable z by z≈Ukz~. Finally, after approximating the nonlinear term as shown in Section 3.2, we obtain the POD-DEIM reduced system in the following form:
(25)M~dz~dt=A~z~+d~+QE~pz~,z~t0=UkTz0,where Q=UkTU¯rPTU¯r−1 and E~pz~=PTEUkz~ as discussed in (22). The next section extends the previous setup for reduced-order modeling for the parametrized miscible flow model.


## 4. Model Reduction for the Parametrized Miscible Flow Model by Using the POD-DEIM Approach withK-Means Clustering

This section appliesK-means clustering with the POD-DEIM reduced-order modeling described in Section 3 on the parametrized miscible flow model given in Section 2. In this work, the numbers of spatial grid points are nx and ny=100, i.e., the dimension of the full-order model is n=nxny. The computation is terminated at the time tf with time stepsize Δt. Moreover, we are interested in an exothermic reaction, which implies that sgnϕ=1. The parameters α, d, Da, and Pe are fixed to be 2, 0.1, 0.01, and 250, respectively. Furthermore, the random noise between 0 and 1 is added to the initial condition for concentration at the collocation points (x^=αPe/2 for all y∈0,Pe) to generate the instability which employs the viscous fingering.


### 4.1.K-Means Clustering

LetIRc⊂ℝ and ILe⊂ℝ be the closed intervals for the values of Rc and Le, respectively. Two parameters Rc,Le are generated randomly in IRc×ILe⊂ℝ2 for 1000 samples and clustered into K centroids using the K-means clustering algorithm as shown in Figure 1.Figure 1
Plot of 1000 random samples for parametersRc,Le∈2,2.5×1,2⊂ℝ2 is clustered into K=20 centroids using the K-means clustering algorithm. The “∗” denotes the centroid of each cluster.TheK-means clustering algorithm for generating centroids is shown in Algorithm 2 as follows. From Algorithm 2, the requirements are number of desire cluster K, number of iteration iter, and the data set which is the set RcLe that contains 1000 training samples of parameters Rc,Le∈2,2.5×1,2⊂ℝ2 in this work. The procedure starts with first randomly picking K initial centroids, p1,p2,⋯,pK. Then, every point in training set measured the distance between each centroid. The points that have a minimum distance to centroid pk are gathered together in the kth cluster. The new centroids will be updated by averaging the whole samples in each cluster. This procedure continues and terminates at the l=iter. As discussed previously, we use K=20 in this work.Algorithm 2:K-means clustering algorithm.
Require: Inputs (i) K∈ℕ, (ii) iter∈ℕ, (iii) RcLe=Rci,LeiT:Rci∈IRc⊂ℝnandLei∈ILe⊂ℝni=11000⊂ℝ2Ensure: Outputs p1,p2,⋯,pK1: Choose randomlyK initial centroids, p1,p2,⋯,pK from RcLe2:D=zeros1000,K+13:forl=1:iterdo4:fori=1:1000do5:forj=1:Kdo6:Di,j=RcLe:,i−pj27:end for8:m,label=minDi,1:K9:Di,K+1=label10:end for11:fork=1:Kdo12:pk⟵meanRcLe:,D:,K+1==k13:end for14:end for


### 4.2. The POD-DEIM Parametrized Reduced Systems

The snapshots for constructing the POD basis are obtained from the full-order models corresponding to all selected centroids from theK-means clustering algorithm. Let zsℓ=ψsℓ,ωsℓ,csℓ,Tsℓ be the numerical solution of zpℓ;ts at time ts from the discretized system of (26) with p=pℓ, which is the ℓth parameter centroid from K-means clustering, i.e., we obtain the snapshots from solving the following 20 full-order systems:
(26)Mpℓdzdt=Apℓz+dpℓ+Epℓ;z,t,ℓ=1,⋯,20.The snapshotsψsℓ,ωsℓ, csℓ, and Tsℓ are the numerical solutions of (26) that approximate ψpℓ;ts, ωpℓ;ts, cpℓ;ts, and Tpℓ;ts, respectively, for s=1,⋯,ns and ℓ=1,⋯K.LetC=c11,c21,⋯,cns1∣c12,c22,⋯,cns2∣⋯∣c1K,c2K,⋯,cnsK∈ℝn×Kns be the snapshot matrix for concentration. In our numerical experiment, K=20 and ns=1500. To obtain the POD basis, the singular value decomposition (SVD) of C is computed, i.e., C=UcΣcVcT. Then, the POD basis of dimension k≪rankC is the first k columns of the left singular matrix Uc of C denoted by Uck∈ℝn×k. Similarly, let UTk, Uψk, and Uωk be the POD basis of dimension k for temperature, stream function, and vorticity, respectively. The POD reduced system is constructed by applying the Galerkin projection on (26) with the basis Uk=diagUψk,Uωk,Uck,UTk∈ℝ4n×4k. The state variable is approximated by zt≈Ukz~t, where z~t=ψ~t,ω~t,c~t,T~t∈ℝ4k is the state variable of the POD reduced-order system, i.e., ct≈Uckc~t, Tt≈UTkT~t, ψt≈Uψkψ~t, and ωt≈Uωkω~t. The resulting POD reduced system is in the form of
(27)M~pdz~dt=A~pz~+d~p+E~p;z~,t,where M~p=UkTMpUk,A~p=UkTApUk∈ℝ4k×4k, and d~p=UkTdp∈ℝ4k are constant matrices for a given parameter vector p, and E~p;z~,t=UkTEp;Ukz~,t is the projected nonlinear term. As described in Section 3, we can obtain a POD-DEIM reduced system from the POD reduced system (27) as shown below:
(28)M~pdz~dt=A~pz~+d~p+E^p;z~,t,where E^p;z~,t is the nonlinear term obtained from using the DEIM approximation for each nonlinear term in E~ in (27). Note that E^p;z~,t is not directly constructed as explained in Section 3, since it consists of many nonlinear functions Nψt,ct,Nψt,Tt,Fψt,ct,Fψt,Tt, and fct given in (13) that depends on different state variables. Therefore, we separately apply this approximation to each nonlinear term. For example, let us first construct the nonlinear snapshot matrix
(29)Nc=N11,N21,⋯,Nns1∣N12,N22,⋯,Nns2∣⋯∣N1K,N2K,⋯,NnsK∈ℝn×Knsfor Nψt,ct, where Nsℓ=Nψsℓ,csℓ. The DEIM is then applied on the nonlinear term Nψt,ct. Other nonlinear terms are done in the same manner. Details on the construction of each nonlinear term by the DEIM can be found in Appendix A which follows from [20]. Finally, the construction of E^ is just the combination of the DEIM results in each nonlinear term. The numerical simulations based on the techniques described in this section are illustrated and analyzed in Section 5.


## 4.1.K-Means Clustering

LetIRc⊂ℝ and ILe⊂ℝ be the closed intervals for the values of Rc and Le, respectively. Two parameters Rc,Le are generated randomly in IRc×ILe⊂ℝ2 for 1000 samples and clustered into K centroids using the K-means clustering algorithm as shown in Figure 1.Figure 1
Plot of 1000 random samples for parametersRc,Le∈2,2.5×1,2⊂ℝ2 is clustered into K=20 centroids using the K-means clustering algorithm. The “∗” denotes the centroid of each cluster.TheK-means clustering algorithm for generating centroids is shown in Algorithm 2 as follows. From Algorithm 2, the requirements are number of desire cluster K, number of iteration iter, and the data set which is the set RcLe that contains 1000 training samples of parameters Rc,Le∈2,2.5×1,2⊂ℝ2 in this work. The procedure starts with first randomly picking K initial centroids, p1,p2,⋯,pK. Then, every point in training set measured the distance between each centroid. The points that have a minimum distance to centroid pk are gathered together in the kth cluster. The new centroids will be updated by averaging the whole samples in each cluster. This procedure continues and terminates at the l=iter. As discussed previously, we use K=20 in this work.Algorithm 2:K-means clustering algorithm.
Require: Inputs (i) K∈ℕ, (ii) iter∈ℕ, (iii) RcLe=Rci,LeiT:Rci∈IRc⊂ℝnandLei∈ILe⊂ℝni=11000⊂ℝ2Ensure: Outputs p1,p2,⋯,pK1: Choose randomlyK initial centroids, p1,p2,⋯,pK from RcLe2:D=zeros1000,K+13:forl=1:iterdo4:fori=1:1000do5:forj=1:Kdo6:Di,j=RcLe:,i−pj27:end for8:m,label=minDi,1:K9:Di,K+1=label10:end for11:fork=1:Kdo12:pk⟵meanRcLe:,D:,K+1==k13:end for14:end for


## 4.2. The POD-DEIM Parametrized Reduced Systems

The snapshots for constructing the POD basis are obtained from the full-order models corresponding to all selected centroids from theK-means clustering algorithm. Let zsℓ=ψsℓ,ωsℓ,csℓ,Tsℓ be the numerical solution of zpℓ;ts at time ts from the discretized system of (26) with p=pℓ, which is the ℓth parameter centroid from K-means clustering, i.e., we obtain the snapshots from solving the following 20 full-order systems:
(26)Mpℓdzdt=Apℓz+dpℓ+Epℓ;z,t,ℓ=1,⋯,20.The snapshotsψsℓ,ωsℓ, csℓ, and Tsℓ are the numerical solutions of (26) that approximate ψpℓ;ts, ωpℓ;ts, cpℓ;ts, and Tpℓ;ts, respectively, for s=1,⋯,ns and ℓ=1,⋯K.LetC=c11,c21,⋯,cns1∣c12,c22,⋯,cns2∣⋯∣c1K,c2K,⋯,cnsK∈ℝn×Kns be the snapshot matrix for concentration. In our numerical experiment, K=20 and ns=1500. To obtain the POD basis, the singular value decomposition (SVD) of C is computed, i.e., C=UcΣcVcT. Then, the POD basis of dimension k≪rankC is the first k columns of the left singular matrix Uc of C denoted by Uck∈ℝn×k. Similarly, let UTk, Uψk, and Uωk be the POD basis of dimension k for temperature, stream function, and vorticity, respectively. The POD reduced system is constructed by applying the Galerkin projection on (26) with the basis Uk=diagUψk,Uωk,Uck,UTk∈ℝ4n×4k. The state variable is approximated by zt≈Ukz~t, where z~t=ψ~t,ω~t,c~t,T~t∈ℝ4k is the state variable of the POD reduced-order system, i.e., ct≈Uckc~t, Tt≈UTkT~t, ψt≈Uψkψ~t, and ωt≈Uωkω~t. The resulting POD reduced system is in the form of
(27)M~pdz~dt=A~pz~+d~p+E~p;z~,t,where M~p=UkTMpUk,A~p=UkTApUk∈ℝ4k×4k, and d~p=UkTdp∈ℝ4k are constant matrices for a given parameter vector p, and E~p;z~,t=UkTEp;Ukz~,t is the projected nonlinear term. As described in Section 3, we can obtain a POD-DEIM reduced system from the POD reduced system (27) as shown below:
(28)M~pdz~dt=A~pz~+d~p+E^p;z~,t,where E^p;z~,t is the nonlinear term obtained from using the DEIM approximation for each nonlinear term in E~ in (27). Note that E^p;z~,t is not directly constructed as explained in Section 3, since it consists of many nonlinear functions Nψt,ct,Nψt,Tt,Fψt,ct,Fψt,Tt, and fct given in (13) that depends on different state variables. Therefore, we separately apply this approximation to each nonlinear term. For example, let us first construct the nonlinear snapshot matrix
(29)Nc=N11,N21,⋯,Nns1∣N12,N22,⋯,Nns2∣⋯∣N1K,N2K,⋯,NnsK∈ℝn×Knsfor Nψt,ct, where Nsℓ=Nψsℓ,csℓ. The DEIM is then applied on the nonlinear term Nψt,ct. Other nonlinear terms are done in the same manner. Details on the construction of each nonlinear term by the DEIM can be found in Appendix A which follows from [20]. Finally, the construction of E^ is just the combination of the DEIM results in each nonlinear term. The numerical simulations based on the techniques described in this section are illustrated and analyzed in Section 5.


## 5. Numerical Experiments

In this section, there are two main numerical tests that consider accuracy and efficiency in terms of reducing the simulation time of the model reduction framework presented in this paper. In Section5.1, the numerical solutions of the proposed method and the original full-order system are compared for different parameter values that are not in the training set. In Section 5.2, the overall computational time reduction and average errors are considered for the solutions from the proposed method when using various randomly selected parameter values in the interested domain. The relevant MATLAB code contributed in this work is available at https://bit.ly/3airQ2u.In this work, the number of spatial grid points isnx=150 and ny=100, i.e., the dimension of the full-order model is n=nxny=15000. The computation is terminated at the time tf=150 with Δt=0.1 (1500 time steps). The system temperature is set to be constant (isothermal case), i.e., it is equivalent to set RT=0.To investigate the possibility of reducing the dimension, we first consider the singular values of both solution snapshots and nonlinear snapshots as depicted in Figure2. Note that this work illustrates only the approximations of concentration. Other state variables have similar accuracy results since they are all coupled in the differential equations.Figure 2
Plots of the first 100 singular values for solution snapshots and nonlinear snapshots.


### 5.1. Numerical Test I

In this section, the POD and POD-DEIM reduced-order model are compared to the full-order model for different given parametersRc,Le∈2,2.5×1,2⊂ℝ2. To demonstrate the efficiency of the reduced-order model, the error at a certain time ts and the global error and its average are defined in terms of the Frobenius norm as follows:
(30)Ets=cts−UckctsF,E=∑s=1nscts−UckctsF,E^=Ens,where cts and c~ts are the concentrations obtained from the full-order and reduced-order models, respectively, at the time step s or the time ts, and ns is the number of snapshots. In this experiment, ns=1500. The approximation c~ts can be the solution from the POD or POD-DEIM reduced-order model. Some cases of parameters Rc,Le∈2,2.5×1,2 are chosen arbitrarily, and the corresponding numerical results are shown below. In particular, we will consider the cases of Rc,Le=2,1,2.5,2, and 2.4,1.1, in Sections 5.1.1, 5.1.2, and 5.1.3, respectively.In Sections5.1.1 and 5.1.2, the parameters are chosen from different corners of our interested 2-dimensional parameter domain, i.e, Rc,Le=2,1 and Rc,Le=2.5,2. These parameters are not in the set of training parameters. As shown in Figures 3 and 4, the resulting POD-DEIM approximation with dimension POD=40/DEIM=60 of concentration for both cases is compared to the full-order model at the time t=100 (s=1000) and the final time t=150 (s=1500). Since these solutions are visually indistinguishable, the absolute errors at all grid points are also provided in Figures 3 and 4. These errors are shown to be in the order of O10−3 or O10−2. More details on the error at specific time t=100, t=150, and average global error are given in Tables 1 and 2, respectively, for Sections 5.1.1 and 5.1.2. As expected, the increasing number of POD basis k provides more accuracy as the error seems to be decreasing. However, this can increase the simulation time as shown in Table 3, which will be discussed in more detail later in Section 5.2.Figure 3
Comparison of approximate concentration for fixed parametersRc,Le=2,1 compared to the full-order model. The results are illustrated at t=100 or s=1000 (top row) and t=150 or s=1500 (bottom row).Figure 4
Comparison of approximate concentration for fixed parametersRc,Le=2.5,2 compared to the full-order model. The results are illustrated at t=100 or s=1000 (top row) and t=150 or s=1500 (bottom row).Table 1
Error at specific time and global error.
Rc,Le=2,1krE100E150EE^FOM (15,000)—————POD40—0.02230.04890.82835.5220×10−4POD60—0.00210.00510.08885.9200×10−5POD-DEIM40400.09730.16503.46072.3071×10−3POD-DEIM40600.03340.06511.23828.2547×10−4Table 2
Error at specific time and global error.
Rc,Le=2.5,2krE100E150EE^FOM (15,000)—————POD40—0.07950.31742.88591.9239×10−3POD60—0.00830.10920.57323.8213×10−4POD-DEIM40400.24910.51009.54236.3621×10−3POD-DEIM40600.03950.35542.36111.5741×10−3Table 3
Average CPU time usage and average global error.
krAvg. CPU timeAvg.EAvg.E^FOM (15,000)——112.8980—POD40—49.28210.79335.2887×10−4POD60—75.53420.13498.9933×10−4POD-DEIM40401.00092.70021.8001×10−4POD-DEIM40601.31270.90246.0160×10−4Similarly, in Section5.1.3, all results are demonstrated in the same manner as seen in Figure 5 and Table 4. The difference is the parameters that are chosen inside the parameter domain, i.e., Rc,Le=2.4,1.1. Table 4 shows that the error at a specific time and the average global error are smaller compared to the previous cases.Figure 5
Comparison of approximate concentration for fixed parametersRc,Le=2.4,1.1 compared to the full-order model. The results are illustrated at t=100 or s=1000 (top row) and t=150 or s=1500 (bottom row).Table 4
Error at specific time and global error.
Rc,Le=2.4,1.1krE100E150EE^FOM (15,000)—————POD40—0.03290.05171.06917.1273×10−4POD60—0.00460.00780.15731.0487×10−4POD-DEIM40400.08520.15302.97591.9839×10−3POD-DEIM40600.02400.03260.81895.4593×10−4


#### 5.1.1.Rc,Le=2,1

#### 5.1.2.Rc,Le=2.5,2

#### 5.1.3.Rc,Le=2.4,1.1

### 5.2. Numerical Test II

In this numerical test, we randomly generate 10 distinct parametersRc,Le in our interested parameter domain 2,2.5×1,2. The solutions are obtained from solving the reduced-order model corresponding to those parameters. CPU time and global error are recorded and averaged as shown in Table 3.From Table3, applying the POD method can make the simulation faster, but it is only approximately twice or three times speedup due to the inefficiency in computing the nonlinear term. After the DEIM is applied, the average CPU time speedups by roughly the order of O100 at equivalent accuracy. In particular, the POD method with k=40 uses approximately 1.3 hours (49,2821 seconds), while the POD-DEIM method with k=40, r=60 uses roughly only 1.3 seconds. In this case, the average error for both POD and POD-DEIM approach is approximately the same. As shown in the previous numerical test, an increasing dimension can improve the accuracy with a trade-off in computational time.The results in this section suggest that the POD-DEIM reduced-order model with the basis constructed from theK-means approach is reliable in approximating the solution for any given parameter in the interested parameter domain 2,2.5×1,2, even though they are not in the training set.


## 5.1. Numerical Test I

In this section, the POD and POD-DEIM reduced-order model are compared to the full-order model for different given parametersRc,Le∈2,2.5×1,2⊂ℝ2. To demonstrate the efficiency of the reduced-order model, the error at a certain time ts and the global error and its average are defined in terms of the Frobenius norm as follows:
(30)Ets=cts−UckctsF,E=∑s=1nscts−UckctsF,E^=Ens,where cts and c~ts are the concentrations obtained from the full-order and reduced-order models, respectively, at the time step s or the time ts, and ns is the number of snapshots. In this experiment, ns=1500. The approximation c~ts can be the solution from the POD or POD-DEIM reduced-order model. Some cases of parameters Rc,Le∈2,2.5×1,2 are chosen arbitrarily, and the corresponding numerical results are shown below. In particular, we will consider the cases of Rc,Le=2,1,2.5,2, and 2.4,1.1, in Sections 5.1.1, 5.1.2, and 5.1.3, respectively.In Sections5.1.1 and 5.1.2, the parameters are chosen from different corners of our interested 2-dimensional parameter domain, i.e, Rc,Le=2,1 and Rc,Le=2.5,2. These parameters are not in the set of training parameters. As shown in Figures 3 and 4, the resulting POD-DEIM approximation with dimension POD=40/DEIM=60 of concentration for both cases is compared to the full-order model at the time t=100 (s=1000) and the final time t=150 (s=1500). Since these solutions are visually indistinguishable, the absolute errors at all grid points are also provided in Figures 3 and 4. These errors are shown to be in the order of O10−3 or O10−2. More details on the error at specific time t=100, t=150, and average global error are given in Tables 1 and 2, respectively, for Sections 5.1.1 and 5.1.2. As expected, the increasing number of POD basis k provides more accuracy as the error seems to be decreasing. However, this can increase the simulation time as shown in Table 3, which will be discussed in more detail later in Section 5.2.Figure 3
Comparison of approximate concentration for fixed parametersRc,Le=2,1 compared to the full-order model. The results are illustrated at t=100 or s=1000 (top row) and t=150 or s=1500 (bottom row).Figure 4
Comparison of approximate concentration for fixed parametersRc,Le=2.5,2 compared to the full-order model. The results are illustrated at t=100 or s=1000 (top row) and t=150 or s=1500 (bottom row).Table 1
Error at specific time and global error.
Rc,Le=2,1krE100E150EE^FOM (15,000)—————POD40—0.02230.04890.82835.5220×10−4POD60—0.00210.00510.08885.9200×10−5POD-DEIM40400.09730.16503.46072.3071×10−3POD-DEIM40600.03340.06511.23828.2547×10−4Table 2
Error at specific time and global error.
Rc,Le=2.5,2krE100E150EE^FOM (15,000)—————POD40—0.07950.31742.88591.9239×10−3POD60—0.00830.10920.57323.8213×10−4POD-DEIM40400.24910.51009.54236.3621×10−3POD-DEIM40600.03950.35542.36111.5741×10−3Table 3
Average CPU time usage and average global error.
krAvg. CPU timeAvg.EAvg.E^FOM (15,000)——112.8980—POD40—49.28210.79335.2887×10−4POD60—75.53420.13498.9933×10−4POD-DEIM40401.00092.70021.8001×10−4POD-DEIM40601.31270.90246.0160×10−4Similarly, in Section5.1.3, all results are demonstrated in the same manner as seen in Figure 5 and Table 4. The difference is the parameters that are chosen inside the parameter domain, i.e., Rc,Le=2.4,1.1. Table 4 shows that the error at a specific time and the average global error are smaller compared to the previous cases.Figure 5
Comparison of approximate concentration for fixed parametersRc,Le=2.4,1.1 compared to the full-order model. The results are illustrated at t=100 or s=1000 (top row) and t=150 or s=1500 (bottom row).Table 4
Error at specific time and global error.
Rc,Le=2.4,1.1krE100E150EE^FOM (15,000)—————POD40—0.03290.05171.06917.1273×10−4POD60—0.00460.00780.15731.0487×10−4POD-DEIM40400.08520.15302.97591.9839×10−3POD-DEIM40600.02400.03260.81895.4593×10−4


### 5.1.1.Rc,Le=2,1

### 5.1.2.Rc,Le=2.5,2

### 5.1.3.Rc,Le=2.4,1.1

## 5.1.1.Rc,Le=2,1

## 5.1.2.Rc,Le=2.5,2

## 5.1.3.Rc,Le=2.4,1.1

## 5.2. Numerical Test II

In this numerical test, we randomly generate 10 distinct parametersRc,Le in our interested parameter domain 2,2.5×1,2. The solutions are obtained from solving the reduced-order model corresponding to those parameters. CPU time and global error are recorded and averaged as shown in Table 3.From Table3, applying the POD method can make the simulation faster, but it is only approximately twice or three times speedup due to the inefficiency in computing the nonlinear term. After the DEIM is applied, the average CPU time speedups by roughly the order of O100 at equivalent accuracy. In particular, the POD method with k=40 uses approximately 1.3 hours (49,2821 seconds), while the POD-DEIM method with k=40, r=60 uses roughly only 1.3 seconds. In this case, the average error for both POD and POD-DEIM approach is approximately the same. As shown in the previous numerical test, an increasing dimension can improve the accuracy with a trade-off in computational time.The results in this section suggest that the POD-DEIM reduced-order model with the basis constructed from theK-means approach is reliable in approximating the solution for any given parameter in the interested parameter domain 2,2.5×1,2, even though they are not in the training set.


## 6. Conclusions

This work provides a model order reduction technique for simulating the parametrized viscous fingering model in a horizontal flow through a porous media domain by using several techniques. TheK-means clustering algorithm is first used to cluster parameters equipped in the model. The snapshots are collected from solving the solution from each centroid. Then, the POD is applied in conjunction with the DEIM to construct a reduced-order model. The resulting reduced-order model can be calculated entirely in the reduced dimension without dependence on the original dimension for various parameter values that are not used in the training set obtained from the K-means algorithm. The approximation admits faster evaluation while trading a small amount of accuracy. The error estimates of the proposed approach for this nonlinear miscible flow model can be obtained by incorporating the accuracy of K-means clustering with the error analysis from [23] and from [24] for, respectively, a posteriori error and a priori error estimates. These theoretical results are beyond the scope of this work and left to be considered in the future.



---

*Source: 3904606-2020-09-04.xml*