The modeling of multiphase flow in a pipe presents a significant challenge for high-resolution computational fluid dynamics (CFD) models due to the high aspect ratio (length over diameter) of the domain. In subsea applications, the pipe length can be several hundreds of meters vs a pipe diameter of just a few inches. Approximating CFD models in a low-dimensional space, reduced-order models have been shown to produce accurate results with a speed-up of orders of magnitude. In this paper, we present a new AI-based non-intrusive reduced-order model within a domain decomposition framework (AI-DDNIROM), which is capable of making predictions for domains significantly larger than the domain used in training. This is achieved by (i) using a domain decomposition approach; (ii) using dimensionality reduction to obtain a low-dimensional space in which to approximate the CFD model; (iii) training a neural network to make predictions for a single subdomain; and (iv) using an iteration-by-subdomain technique to converge the solution over the whole domain. To find the low-dimensional space, we compare Proper Orthogonal Decomposition with several types of autoencoder networks, known for their ability to compress information accurately and compactly. The comparison is assessed with two advection-dominated problems: flow past a cylinder and slug flow in a pipe. To make predictions in time, we exploit an adversarial network, which aims to learn the distribution of the training data, in addition to learning the mapping between particular inputs and outputs. This type of network has shown the potential to produce visually realistic outputs. The whole framework is applied to multiphase slug flow in a horizontal pipe for which an AI-DDNIROM is trained on high-fidelity CFD simulations of a pipe of length 10 m with an aspect ratio of 13:1 and tested by simulating the flow for a pipe of length 98 m with an aspect ratio of almost 130:1. Inspection of the predicted liquid volume fractions shows a good match with the high fidelity model as shown in the results. Statistics of the flows obtained from the CFD simulations are compared to those of the AI-DDNIROM predictions to demonstrate the accuracy of our approach.

## I. INTRODUCTION

Non-intrusive reduced-order modeling (NIROM) has been the subject of intense research activity over the last five years, largely due to the advances made in machine learning and the re-application of these techniques to reduced-order modeling. This paper takes a step toward demonstrating how non-intrusive reduced-order models can generalize by training the model on one domain and deploying it on a much larger domain. The method outlined here could be extremely useful for the energy industry, for example, in which pipelines of the order of kilometers in lengths and inches in diameter are used for the subsea transportation of fluids. With such high aspect ratios, these pipes are too long to be modeled by high-resolution computational fluid dynamics (CFD) models alone. In this paper, we propose a non-intrusive reduced-order model based on autoencoders (for dimensionality reduction), an adversarial network (for prediction), and a domain decomposition approach. For dimensionality reduction, we investigate the performance of several autoencoders and Proper Orthogonal Decomposition (POD) using two test cases (flow past a cylinder and multiphase slug flow in a horizontal pipe). With the method that performs best (the convolutional autoencoder), we go on to demonstrate the NIROM approach on multiphase slug flow in a horizontal pipe, training the networks on CFD data from a 10 m pipe with an aspect ratio of 13:1 and making predictions for the flow within a 98 m pipe with an aspect ratio of almost 130:1. In the following paragraphs, we give some background on reduced-order modeling (ROM) and NIROM; on dimensionality reduction methods and the use of autoencoders; on prediction; domain decomposition methods and ROM; and finally on multiphase flow. The final two paragraphs summarize the main contributions of the paper and describe the layout of the rest of the paper.

The aim of reduced-order modeling^{1} (ROM) is to obtain a low-dimensional approximation of a computationally expensive high-dimensional system of discretized equations, henceforth referred to as the high-fidelity model (HFM). To be of benefit, the low-dimensional model should be accurate enough for its intended purpose and orders of magnitude faster to solve than the HFM. Known as projection-based ROM,^{2} one common strategy for constructing reduced-order models is to use a Galerkin projection of the HFM onto a low-dimensional subspace. However, in this article we focus on an alternative method, NIROM, which, unlike projection-based ROM, does not require access to or modification of the source code of the HFM. It requires only the results of the HFM with which it constructs a low-dimensional approximation to the HFM in two stages: the offline stage and the online stage. During the offline stage, solutions from the HFM are generated (known as snapshots); a set of basis functions that span the low-dimensional or reduced space are obtained by a dimensionality reduction method; and finally, the evolution of the HFM in the reduced space is approximated in some manner. The latter step can be done in several ways, but here, as we focus on AI-based NIROM, we use a neural network. A profusion of terms exists for this type of non-intrusive modeling, including POD with interpolation;^{3} NIROM;^{4,5} POD surrogate modeling;^{5,6} system or model identification;^{7,8} Galerkin-free;^{9} data-driven reduced-order modeling;^{10–12} Deep Learning ROM;^{13} and digital twins.^{14–17} In addition to making predictions, digital twins assimilate data from observations to improve the accuracy of the prediction.

To find the low-dimensional subspace in which to approximate the HFM, many of these non-intrusive approaches rely on Proper Orthogonal Decomposition (POD),^{18} which is based on Singular Value Decomposition. Also known as Principal Component Analysis, POD finds the optimal linear subspace (with a given dimension) that can represent the space spanned by the snapshots and prioritizes the modes according to those that exhibit the most variance. Whilst POD works well in many situations, for advection-dominated flows with their slow decay of singular values or large Kolmogorov *N*–width,^{19} approximations based on POD can be poor^{20–22} and researchers are turning increasingly to autoencoders.^{23} Although adding to the offline cost, these networks seek a low-dimensional nonlinear subspace, which can be more accurate and efficient than a linear subspace for approximating the HFM.

Convolutional networks are particularly good at analyzing and classifying images (on structured grids)^{24,25} with the ability to pick out features and patterns wherever their location (translational invariance), and these methods are applicable directly to the dimensionality reduction of CFD solutions on structured grids through the use of convolutional autoencoders (CAEs). Methods that apply convolutional networks to data on unstructured meshes do exist (based on space-filling curves;^{26} graph convolutional networks;^{27,28} and a method that introduces spatially varying kernels^{29}) but are in their infancy, so most researchers either solve the high-resolution problem on structured grids directly or interpolate from the high-fidelity model snapshots to a structured grid before applying the convolutional layers. The latter approach is adopted here.

Perhaps the first use of an autoencoder for dimensionality reduction within a ROM framework was applied to reconstruct flow fields in the near-wall region of channel flow based on information at the wall,^{30} whilst the first use of a convolutional autoencoder came 16 years later and was applied to Burgers Equation, advecting vortices and lid-driven cavity flow.^{31} In the few years since 2018, many papers have appeared, in which convolutional autoencoders have been applied to sloshing waves, colliding bodies of fluid and smoke convection;^{32} flow past a cylinder;^{33–35} the Sod shock test and transient wake of a ship;^{36} air pollution in an urban environment;^{37–39} parametrized time-dependent problems;^{40} natural convection problems in porous media;^{41} the inviscid shallow water equations;^{42} supercritical flow around an airfoil;^{43} cardiac electrophysiology;^{44} multiphase flow examples;^{45} the Kuramoto–Sivashinsky equation;^{46} the parametrized 2D heat equation;^{47} and a collapsing water column.^{48} Of these papers, those which compare autoencoder networks with POD generally conclude that autoencoders can outperform POD,^{31,33} especially when small numbers of reduced variables are used.^{41–44} However, when large enough numbers of POD basis functions are retained, POD can yield good results, sometimes outperforming the autoencoders.

A recent dimensionality reduction method that combines POD/SVD and an autoencoder (SVD-AE) has been introduced independently by a number of researchers and demonstrated on: vortex-induced vibrations of a flexible offshore riser at a high Reynolds number^{49} (described as hybrid ROM); the generalized eigenvalue problems associated with neutron diffusion^{50} (described as an SVD autoencoder); Marsigli flow^{51} (described as nonlinear POD); and cardiac electrophysiology^{52} (described as POD-enhanced deep learning ROM). This method has at least three advantages: (i) by training the autoencoder with POD coefficients, it is of no consequence whether the snapshots are associated with a structured or unstructured mesh; (ii) an initial reduction of the number of variables by applying POD means that the autoencoder will have fewer trainable parameters and therefore be easier to train; and (iii) autoencoders in general can find the minimum number of latent variables needed in the reduced representation. For example, the solution of flow past a cylinder evolves on a one-dimensional manifold parametrized by time; therefore, only one latent variable is needed to capture the physics of this solution.^{26,42,44}

The Adversarial Autoencoder^{53} (AAE) is a generative autoencoder sharing similarities with the variational autoencoder (VAE) and the generative adversarial network (GAN). In addition to an encoder and decoder, the AAE has a discriminator network linked to its bottleneck layer. The purpose of the discriminator and associated adversarial training is to make the posterior distribution of the latent representation close to an arbitrary prior distribution thereby reducing the likelihood that the latent space will have “gaps.” Therefore, any set of latent variables should be associated, through the decoder, with a visually realistic output. Not many examples exist of using an AAE for dimensionality reduction in fluid dynamics problems; however, it has been applied to model air pollution in an urban environment.^{38,39} In this work, we compare POD, CAE, AAE, and the SVD-AE on flow past a cylinder and multiphase flow in a pipe, to assess their suitability as dimension reduction methods.

Once the low-dimensional space has been found, the snapshots are projected onto this space, and the resulting reduced variables (either POD coefficients or latent variables of an autoencoder) can be used to train a neural network, which attempts to learn the evolution of the reduced variables in time (and/or their dependence on a set of parameters). From the references in this paper alone, many examples exist of feed-forward and recurrent neural networks having been used for the purpose of learning the evolution of time series data, for example, by Multi-layer perceptrons,^{12,13,40,41,43,54–60} Gaussian Process Regression,^{11,45,61–63} and Long-Short Term Memory networks.^{31,32,34,35,38,51,64} When using these types of neural networks to predict in time, if the reduced variables stray outside of the range of values encountered during training, the neural network can produce unphysical, divergent results.^{39,51,52,64,65} To combat this, a number of methods have been proposed. Physics-informed neural networks^{55} aim to constrain the predictions of the neural network to satisfy physical laws, such as conservation of mass or momentum.^{59,60} A method introduced by Refs. 56 and 57 aims to learn the mapping from the reduced variables at a particular time level to their time derivative, rather than the reduced values themselves at a future time level. This enables the use of variable time steps when needed, to control the accuracy of the solution in time. A third way of tackling this issue, which is explored in this paper, is to use adversarial networks, renowned for their ability to give realistic predictions.

Adversarial networks, such as the GAN and the AAE, aim to learn a distribution to which the training data could belong, in addition to a mapping between solutions at successive time levels. GANs and AAEs are similar in that they both use a discriminator network and deploy adversarial training, and both require some modification so that they can make predictions in time. The aim of these networks is to generate images (or in this case, reduced variables associated with fluid flows) that are as realistic as possible. To date, there are not many examples of the use of GANs or AAEs for prediction in CFD modeling. Two exceptions are Ref. 66, which combines a VAE and GAN to model flow past a cylinder and the collapse of a water dam and Ref. 67, which uses a GAN to predict the reduced variables of an epidemiological model which modeled the spread of a virus through a small, idealized town. This particular model performed well when compared with an LSTM.^{68} Conditional GANs (CGAN) have similar properties to the GAN and AAE, and they have been used successfully to model forward and inverse problems for coupled hydro-mechanical processes in heterogeneous porous media;^{69} a flooding event in Hokkaido, Japan, after the 1993 earthquake;^{70} and a flooding event in Denmark.^{71} However, the closeness of the CGAN's distribution to that of the training data is compromised by the “condition” or constraint. GANs are known to be difficult to train, so, in this paper, we use an Adversarial Autoencoder, albeit modified, so that it can predict the evolution of the reduced variables in time.

Combining domain decomposition techniques with ROM has been done by a number of researchers. An early example^{72} presents a method for projection-based ROMs in which the POD basis functions are restricted to the nodes of each subdomain of the partitioned domain. A similar approach has also been developed for non-intrusive ROMs,^{61} which was later extended to partition the domain by minimizing communication between subdomains,^{62} effectively isolating as much as possible, the physical complexities between subdomains. As the domain of our main test case (multiphase flow in a pipe) is long and thin with a similar amount of resolution and complexity of behavior occurring in partitions of equal length in the axial direction, here, we simply split the domain into subdomains of equal length in the axial direction (see Fig. 1). The neural network learns how to predict the solution for a given subdomain, and the solution throughout the entire pipe is built up by using the iteration-by-subdomain approach.^{73} The domain decomposition approach we use has some similarities to the method employed in Ref. 74, which decomposes a domain into patches to make training a neural network more tractable. However, our motivation for using domain decomposition is to make predictions for domains that are significantly larger than those used in the training process. When modeling a pipe that is longer than the pipe used to generate the training data, it is likely that the simulation will need to be run for longer than the original model as the fluid will take longer to reach the end of the pipe. This means that boundary conditions for the longer pipe must be generated somehow, rather than relying on using boundary conditions from the original model. Generating suitable boundary conditions for turbulent CFD problems is, in general, an open area of research. Often used are incoming synthetic-eddy methods,^{75} which attempt to match specified mean flows and Reynolds stresses at the inlet. Recently, researchers have explored using GANs to generate boundary conditions with success.^{76,77} We present three methods of generating boundary conditions for our particular application and also discuss alternative methods in Conclusions and Further Work.

The test case of multiphase flow in a pipe is particularly challenging due to the difficulties such as the space-time evolution of multiphase flow patterns (stratified, bubbly, slug, annular), the turbulent phase-to-phase interactions, the drag, inertia, and wake effects that arise for the HFM from the high aspect ratio (length to diameter) of the domain of a typical pipe. Many address this by developing one dimensional (flow regime-dependent or -independent) models for long pipes.^{78–80} Nevertheless, such models contain some uncertainties as they rely on several closure or empirical expressions^{81} under the limited experimental data^{82} in describing, for example, the 3D space-time variations of interfacial frictional forces with phase distributions (the bubble/drop entrainment, the bubble-induced turbulence, and the phase interfacial interactions), depending on the flow pattern, flow direction, and pipe physical properties (inclination, diameter, and length). Significant progress has been made in 3D modeling^{83} by using direct numerical simulations^{84} (DNS) and front-tracking methods.^{85} To generate the solutions of the HFM, we employ a method based on Large Eddy Simulation, which advects a volume fraction field^{86} and uses mesh adaptivity to have high resolution where most needed. Although compromising on resolving features on the smaller temporal and spatial scales, this approach is computationally more feasible than DNS and has the advantage of being conservative, unlike front-tracking methods.

In this paper, we propose a non-intrusive reduced-order model (AI-DDNIROM) capable of making predictions for a domain to which it has not been exposed during training. Several autoencoders are explored for the dimensionality reduction stage, as there is evidence that they are more efficient than POD for advection-dominated problems such as those tackled here. The dimensionality reduction methods are applied to 2D flow past a cylinder and 3D multiphase slug flow in a horizontal pipe. For the prediction stage, an adversarial network is chosen (based on a modified adversarial autoencoder) as these types of networks are believed to generate latent spaces with no gaps^{53} and thus are likely to produce more realistic results than feed-forward or recurrent neural networks without adversarial layers. A domain decomposition approach is applied, which, with an iteration-by-subdomain technique, enables predictions to be made for multiphase slug flow with a significantly longer pipe than was used when training the networks. The predictions of the adversarial network are taken from a probability distribution learned during training. Any point within the Gaussian distribution of the latent variables should therefore result in a realistic solution. Statistics from the HFM solutions and predictions of the non-intrusive reduced-order models for the original length pipe and the longer pipe are compared. The contributions of this work are: (i) a method, which can make predictions for a domain significantly larger than that used to train the reduced-order models; (ii) the exploitation of an adversarial network to make realistic predictions, and comparing statistics of the reduced-order models with the original CFD model; and (iii) the investigation of a number of methods to generate boundary conditions for the larger domain.

The outline of the remainder of the paper is as follows. Section II describes the methods used in constructing the reduced-order models and the domain decomposition approach, which is exploited in order to be able to make predictions for a longer domain than that used in training. Section III presents the results for the dimensionality reduction methods applied to flow past a cylinder and multiphase flow in a pipe and then shows the predictions of the reduced-order model of multiphase flow in a pipe, for both the original domain and the extended domain. Conclusions are drawn, and future work described in the final section. Details of the hyperparameter optimization process and the network architectures are given in the Appendix.

## II. METHODOLOGY

### A. Non-intrusive reduced-order models

The offline stage of a non-intrusive reduced-order model can be split into three steps: (i) generating the snapshots by solving a set of discretized governing equations (the high-resolution or high-fidelity model); (ii) reducing the dimensionality of the discretized system; and (iii) teaching a neural network to predict the evolution of the snapshots in reduced space. The online stage consists of two steps: (i) predicting values of the reduced variables with the neural network for an unseen state and (ii) mapping back to the physical space of the high-resolution model. In this section, the methods used in this investigation for dimensionality reduction (Sec. II B) and prediction (Sec. II C) are described. The final section (Sec. II D) outlines an approach for making predictions for a larger domain having used a smaller domain to generate the training data.

### B. Dimensionality reduction methods

Described here are four techniques for dimensionality reduction, which are used in this investigation, namely Proper Orthogonal Decomposition, a convolutional autoencoder, an adversarial autoencoder, and a hybrid SVD autoencoder.

#### 1. Proper orthogonal decomposition

Proper Orthogonal Decomposition is a commonly used technique for dimensionality reduction when constructing reduced-order models. POD requires the minimization of the reconstruction error of the projection of a set of solutions (snapshots) onto a number of basis functions, which define a low-dimensional space. In order to minimize the reconstruction error, the basis functions must be chosen as the left singular vectors of the singular value decomposition (SVD) of the matrix of snapshots. Suppose the snapshots matrix is represented by ** S**, whose columns are solutions at different instances in time (i.e., the snapshots) and whose rows correspond to nodal values of solution variables, then

**can be decomposed as**

*S*where the matrix ** U** contains the left singular vectors,

**the right singular vectors, and Σ contains the singular values on its diagonal, zeros elsewhere. If POD is well suited to the problem, many of the singular values will be close to zero and the corresponding columns of**

*V***can be discarded. The POD basis functions to be retained are stored in a matrix denoted by**

*U***. The POD coefficients of a snapshot can be found by pre-multiplying the snapshot by $RT$, and the reconstruction of a snapshot can be found by pre-multiplying the POD coefficients of the snapshot by**

*R***:**

*R*where $uk$ is the *k*th snapshot and $(urecon)k$ is its reconstruction. Hence, the reconstruction error over a set of *N* snapshots ${u1,u2,\u2026,uN}$ can be written as

Often the mean is subtracted from the snapshots before applying singular value decomposition; however, in this study, doing so was found to have little effect. In the first test case, 2D flow past a cylinder, two velocity components are included in the snapshot matrix,

where *u _{i}* and

*v*represent the

_{i}*x*and

*y*components of velocity, respectively, at the

*i*th node;

*k*denotes a particular snapshot; and

*M*is the number of nodes. For the 3D multiphase flow test case, the snapshots comprise velocities and volume fractions, so a single snapshot has the form

where $wik$ and $\alpha ik$ represent the *z* component of velocity and the volume fraction, respectively, at the *i*th node of the *k*th snapshot. In this case, the velocity components are scaled to be in the range $[\u22121,1]$ so that their magnitudes are similar to those of the volume fractions.

#### 2. Convolutional autoencoder

An autoencoder is a particular type of feed-forward network that attempts to learn the identity map.^{87} When used for compression, these networks have a central or bottleneck layer that has fewer neurons than the input and output layers, thereby forcing the autoencoder to learn a compressed representation of the training data. An autoencoder consists of an encoder, which compresses the data to the latent variables of the bottleneck layer, and a decoder, which decompresses or reconstructs the latent variables to an output layer of the same dimension as the input layer. The latent variables span what is referred to as the latent space. The convolutional autoencoder typically uses a series of two types of layers to compress the input data in the encoder: convolutional layers and pooling layers. These layers both apply operations to an input grid resulting in an output grid (or feature map) of reduced size. The inverse operations are then used in succession in a decoder, resulting in a reconstructed grid of the same shape as the input. The encoder-decoder pair can be trained as any other neural network: by passing training data through the network and updating the weights associated with the layers according to a loss function such as the mean square error. If $uk$ represents the *k*th sample in the dataset of *N* samples and $(urecon)k$ represents the corresponding output of the autoencoder, which can be written as

then the mean square error can be expressed as in Eq. (3).

#### 3. Adversarial autoencoder

The adversarial autoencoder^{53} is a recently developed neural network that uses an adversarial strategy to force the latent space to follow a (given) prior distribution ($Pprior$). Its encoder-decoder network is the same as that of a standard autoencoder; however, in addition, the adversarial autoencoder includes a discriminator network, which is trained to distinguish between true samples (from the prior) and fake samples (from the latent space). There are therefore three separate training steps per mini-batch. In the first step, the reconstruction error of the inputs is minimized (as is done in a standard autoencoder). In the second and third steps, the adversarial training takes place. In the second step, the discriminator network is trained on latent variables sampled from the prior distribution with label 1 and latent variables generated by the encoder with label 0. In the third step, the encoder is trained to fool the discriminator, that is, it tries to make the discriminator produce an output of 1 from its generated latent vectors. Note that this is the role of the generator in a GAN and, as such, the encoder (*G*) and discriminator (*D*) play the minimax game described by Eq. (7). This equation is the implicit loss function for the adversarial training:

where *V* is the value function that *G* and *D* play the minimax game over, $z\u223cPprior$ is a sample from the desired distribution, and $u\u223cPdata$ is a sample input grid. There are strong similarities between the adversarial autoencoder, GANs and Variational Autoencoders (VAEs). All three types of networks set out to obtain better generalization than non-adversarial networks by attempting to obtain a smooth latent space with no gaps. Results in Ref. 53 show that the AAE performs better at this task than the VAE on the MNIST digits. Imposing a prior distribution upon the variables of the latent space ensures that any set of latent variables, when passed through the decoder, should have a realistic output.^{53}

#### 4. SVD autoencoder

As the name suggests, the SVD autoencoder makes use of two strategies. Initially, an SVD is applied to the data, resulting in POD coefficients that are subsequently used to train an autoencoder, which applies a second level of compression. Once trained, the latent variables of the SVD autoencoder can be written as

where $fenc$ is the encoder, ** R** represents the POD basis functions, $uk$ is the

*k*th snapshot, and $zk$ are the latent variables. For reconstruction, the inverse of this process is then employed, whereby a trained decoder first decompresses the latent space variables to POD coefficients, after which these POD coefficients are reconstructed to the original space of the high-fidelity model. The reconstruction can be written as

where $fdec$ is the decoder, $fae$ is the autoencoder, and $(urecon)k$ is the reconstruction of the *k*th snapshot. This network could be approximated by adding to the autoencoder a linear layer after the input and before the output, and dispensing with the SVD, however, it has been found that this network is harder to train. Here, we take advantage of the efficiency of the SVD and use this in conjunction with an autoencoder.

### C. Prediction methods

In this study, when predicting, we wish to approximate a set of reduced variables (either POD coefficients or latent variables of an autoencoder) at a future time step. The adversarial autoencoder is re-purposed for this task in an attempt to capitalize on the fact that this network should produce realistic results (providing that the training dataset is representative of the behavior that will be modeled). So that it can predict time series data, three modifications are made to the original adversarial autoencoder network:^{53} namely that (i) the bottleneck layer no longer has fewer variables than the input (to prevent further compression); (ii) the output is the network's approximation of the reduced variables at a future time level; and (iii) the input is the reduced variables at the preceding time level as well as the reduced variables of the neighboring subdomains at the future time (as we adopt a domain decomposition approach, which is described in the next paragraph). The modified adversarial autoencoder is trained by minimizing the error between its output and the predicted variables at the future time level, as well as incorporating the adversarial training strategy described in Sec. II B 3. To avoid confusion, we refer to this network as a predictive adversarial network, because, with different inputs and outputs, it is no longer an autoencoder.

In this study, we adopt a domain decomposition approach to facilitate predicting the solution for larger domains than that used in training (see Sec. II D). Given the aspect ratio of the pipe, we split the domain into subdomains of equal length in the axial direction, see Fig. 1. To train the predictive adversarial network, reduced variables are obtained by interpolating the high-fidelity solutions or snapshots onto a structured grid in each subdomain in turn and compressing the interpolated snapshots from all the subdomains using POD or an autoencoder. The interpolation is linear and achieved by using the finite element basis functions. The predictive adversarial network is taught to predict the reduced variables in a particular subdomain at a future time level given the reduced variables in the neighboring subdomains at the future time level and the reduced variables in the subdomain at the preceding time level. Using training data for all the subdomains and those time levels that are in the training dataset, the predictive adversarial network learns the mapping *f*, written as

where $zik$ represents the reduced variables in subdomain *i* at the future time level *k*; $zik\u22121$ represents the same but at the preceding time level; and $zi\u22121k$ and $zi+1k$ denote the reduced variables at the future time level for the subdomains to the left and right of subdomain *i*. When predicting for one time level, all subdomains are iterated over (the iteration-by-subdomain method) until convergence is reached over the whole domain. This is done by sweeping from left to right (increasing *i*) and then sweeping from right to left (decreasing *i*). During the iteration process, $zik$ of Eq. (10) is being continually updated. As we consider incompressible flows in this study, the solution method has to be implicit in order to allow information to travel throughout the domain within one time level. This sweeping from left to right and back again allows information to pass from the leftmost to the rightmost subdomains and vice versa.

For each new time level, an initial solution is required to start the iteration process [for $zi\u22121k$ and $zi+1k$ in Eq. (10)]. The solution at the previous, converged time level could be used ($zi\xb11k=zi\xb11k\u22121$); however, using linear extrapolation based on two previous time levels showed better convergence

The procedure for sweeping over the subdomains is given in Algorithm 1, in which *f* represents the predictive adversarial network, $Ntime$ is the number of time levels, $Nsweep$ is the number of sweeps carried out over the whole domain, and $Nsub$ is the total number of subdomains. Two of these subdomains are treated as boundary conditions and are fully imposed throughout the duration of the prediction, so at line 18 of Algorithm I, only the subdomains where a solution is sought are iterated over. In this study, a fixed number of sweep iterations were used as this gave good results; however, a convergence criterion could be easily implemented if desired.

1: !! set initial conditions for each subdomain i |

2: $zi0\u2009\u2200i$ |

3: for time level $k=1,2,\u2026,Ntime$ do |

4: !! set boundary conditions |

5: $z1k,\u2009zNsubk$ |

6: !! estimate the solution at the future time level k for all the subdomains |

7: if k > 1 then |

8: for subdomain $i=2,3,\u2026,Nsub\u22121$ do |

9: $zik=zik\u22121+(zik\u22121\u2212zik\u22122)$ |

10: end for |

11: else |

12: for subdomain $i=2,3,\u2026,Nsub\u22121$ do |

13: $zik=zik\u22121$ |

14: end for |

15: end if |

16: !! sweep over subdomains |

17: for sweep iteration $j=1,2,\u2026,Nsweep$ do |

18: for subdomain $i=2,3,\u2026,Nsub\u22122,Nsub\u22121,Nsub\u22122,\u2026,4,3$ do |

19: !! calculate the latent variables of subdomain i at time level k |

20: $zik=f(zi\u22121k,\u2009zik\u22121,\u2009zi+1k)$ |

21: end for |

22: end for |

23: end for |

1: !! set initial conditions for each subdomain i |

2: $zi0\u2009\u2200i$ |

3: for time level $k=1,2,\u2026,Ntime$ do |

4: !! set boundary conditions |

5: $z1k,\u2009zNsubk$ |

6: !! estimate the solution at the future time level k for all the subdomains |

7: if k > 1 then |

8: for subdomain $i=2,3,\u2026,Nsub\u22121$ do |

9: $zik=zik\u22121+(zik\u22121\u2212zik\u22122)$ |

10: end for |

11: else |

12: for subdomain $i=2,3,\u2026,Nsub\u22121$ do |

13: $zik=zik\u22121$ |

14: end for |

15: end if |

16: !! sweep over subdomains |

17: for sweep iteration $j=1,2,\u2026,Nsweep$ do |

18: for subdomain $i=2,3,\u2026,Nsub\u22122,Nsub\u22121,Nsub\u22122,\u2026,4,3$ do |

19: !! calculate the latent variables of subdomain i at time level k |

20: $zik=f(zi\u22121k,\u2009zik\u22121,\u2009zi+1k)$ |

21: end for |

22: end for |

23: end for |

### D. Extending the domain

In this study, we investigate the ability of a non-intrusive reduced-order model in combination with a domain decomposition approach to be able to make predictions for domains larger than that used in the training process. We test this approach on the dataset generated from multiphase flow in a pipe. With sufficient initial conditions and boundary conditions, exactly the same procedure can be used to make predictions for the extended domain as is used to make predictions for the domain used in training. That is, the solution is obtained for a single subdomain, whilst sweeping over all subdomains until convergence is reached (outlined in Sec. II C).

As the length of the pipe of interest (“extended pipe”) is longer than the pipe used in training, initial conditions must be generated throughout the extended pipe. The method used here is to specify initial conditions throughout the extended pipe by repeating initial conditions from the shorter pipe. An alternative would be to find the reduced variables for a steady state (for example, water in the bottom half of the pipe and air in the top half) and use these values in every subdomain in the extended pipe. We choose the former method to reduce the time taken for instabilities and slugs to develop.

For the extended pipe, boundary conditions (effectively the reduced variables in an entire subdomain) can be imposed using the data already available from the HFM. However, as the length of the pipe is longer than the pipe used in training, we wish to make predictions over a longer period over which snapshots were collected from the HFM. In order to obtain boundary conditions for the extended pipe, several methods are explored for the inlet or upstream boundary. Of those investigated, three methods performed better than the others, listed below.

- (i)
Cycling through slug formation: a slug is found in the shorter pipe, and the velocity and volume fraction fields associated with the advection of the slug through a subdomain are looped over in the upstream boundary subdomain.

- (ii)
Perturbed instability: the volume fraction field associated with an instability from the shorter pipe is perturbed with Gaussian noise. This is then imposed on the boundary subdomain. The associated velocity field is used unperturbed.

- (iii)
Original boundaries repeated: solutions from the shorter pipe are cycled through in the boundary subdomain.

At the downstream boundary, reduced variables corresponding to a steady state solution (water in the bottom half of the pipe and air in the top half) were imposed. Specific details of the boundary conditions are given in the results section. These three approaches are somewhat heuristic. As we are using information that the model will not have seen and that does not accurately satisfy the governing equations, we exploit the ability of the predictive adversarial network to produce realistic results, as it should have learnt appropriate spatial and temporal covariance information during training. An alternative method for generating boundary conditions is discussed in the section on conclusions and future work.

## III. RESULTS

### A. Test cases

Two test cases are used to demonstrate the dimensionality reduction methods proposed in this paper. The first is flow past a cylinder in 2D; the second is 3D multiphase flow in a pipe. The second test case is also used to demonstrate the prediction capabilities of the predictive adversarial network for both the domain that was used in training and a domain that is significantly longer that the one used in training. The code used to produce the high-fidelity model results, IC-FERST, has been validated with experimental results for a 2D solitary wave and 3D collapsing water column,^{88} and with the Rayleigh–Taylor benchmark problem in 2D and 3D.^{86} The test cases used in this paper are now described.

#### 1. Flow past a cylinder

The following partial differential equations describe the motion of an incompressible fluid:

where *ρ* is the density (assumed constant), ** u** is the velocity vector, $\tau $ contains the viscous terms associated with an isotropic Newtonian fluid,

*p*represents the non-hydrostatic pressure,

*t*is time, and the gradient operator $\u2207$ is defined as

When solving these equations, a linear triangular element is used with a discontinuous Galerkin discretization for the velocities and a continuous Galerkin representation of the pressure (the P1DG-P1CV element). As well as satisfying the Ladyzhenskaya–Babuška–Brezzi condition, this element type is stable and accurate even on highly distorted elements, such as those which may occur along the interface of the two fluids.^{88} To discretise in time, Crank–Nicolson is used. The resulting velocity solutions will satisfy the discretized continuity equation. As the velocity field fully describes incompressible flow, only the velocity variables are required by the reduced-order models. For more details on how this system of equations is discretized and solved, the reader is referred to Ref. 89. For the flow past a cylinder test case, the domain measures 2.2 m (horizontal axis) by 0.41 m (vertical axis), and the center of the cylinder is located at 0.2 m from the leftmost boundary on the horizontal centerline of the domain. Free slip and no normal flow boundary conditions are applied on the upper and lower walls; no slip is applied on the surface of the cylinder. Zero shear and zero normal stress are applied at the outlet (the right-hand boundary of the domain). In the following results, speeds and velocities are given in meters per second and time is in seconds. A Reynolds number of 3900 was used:

where *U* is the constant inlet velocity, $U=0.039\u2009ms\u22121$, the density has value $\rho =1000\u2009kgm\u22123$, and the diameter of the cylinder is $L=0.1\u2009m$. Thus, the dynamic viscosity is $\mu =10\u22123\u2009kg\u2009m\u22121\u2009s\u22121$. Formed from solutions of this problem, the dataset consists of 2000 snapshots with a time interval of 0.25 s. (An adaptive time step was used to solve the equations; however, the solutions were saved every 0.25 s to generate the snapshots.)

#### 2. Multiphase flow in a pipe

Multiphase slug flow in a horizontal pipe is used as the second test case. We use an interface capturing method, in which we track the interface by solving an advection equation for the volume fraction of the liquid phase. Let *α* be the volume fraction of the liquid (water in this case), which means that the volume fraction of the gas (air) is $(1\u2212\alpha )$. The conservation of mass for incompressible fluids can therefore be written as

where *t* represents time and ** u** represents velocity. Assuming incompressible viscous fluids, conservation of momentum yields the following:

where *p* represents pressure, ** g** is the gravitational acceleration vector $(0,0,9.8)$, and $F\sigma $ is the force representing surface tension. The bulk density and bulk viscosity are defined as

respectively, where $\rho water$ and $\rho air$ are the densities of water and air, respectively, and $\mu water$ and $\mu air$ are the dynamic viscosities of water and air, respectively. Again, the resulting velocity solutions will satisfy the discretized continuity equation. For more details of how the governing equations are discretized and solved, see Ref. 86, including information on the unstructured adaptive meshing process, the adaptive time stepping and compressive advection technique to keep the interface at the boundary of the fluids sharp. The densities of air and water are taken as 1.125 and 1000 kg m^{−3}, respectively, and the viscosities are $1.81\xd710\u22125$ and $9.892\xd710\u22124\u2009kg\u2009m\u22121\u2009s\u22121$, respectively. The modeled pipe has dimensions of 10 m in length and a radius of 0.039 m. Boundary conditions of no normal flow and no slip were weakly enforced on the pipe wall, and any incoming momentum is given a value of zero. The outlet of the pipe has a non-hydrostatic pressure of zero, and again, any incoming velocities are set to zero and incoming volume fraction is taken to be water. Initially, the pipe is filled entirely with water, which flows along the axial direction at a velocity of 4.162 m s^{−1} in the top half of the pipe and 2.082 m s^{−1} in the bottom half. After the first time step, air starts flowing in through the inlet through the top half at a velocity of 4.162 m s^{−1}, a scenario which can lead to the formation of slugs. These values of velocity correspond to superficial velocities of air and water of 2.081 and 1.041 m s^{−1}, respectively. The dataset used for training the reduced-order models consists of solutions at 800 time levels with a fixed time interval of 0.01 s. The reduced-order models use the velocity fields in three directions and the volume fraction field.

### B. Dimensionality reduction

Four methods for dimensionality reduction (or compression) are compared, namely POD, CAE, AAE, and SVD-AE. An extensive hyperparameter optimsation was performed to find the optimal set of values for the hyperparameters of each autoencoder. Details of the hyperparameters that were varied, the ranges over which they were varied, and the optimal values and architectures that were obtained as a result can be found in Tables IV–,VI. Ten POD basis functions were retained for the compression based on POD and ten latent variables were used in the bottleneck layers of the autoencoders. For the SVD-AE, one hundred POD coefficients were retained, which were then compressed to ten latent variables by an autoencoder. The top part (shaded blue) of Fig. 2 shows a schematic diagram of how the networks used for dimensionality reduction are trained for the flow past a cylinder test case.

#### 1. Flow past a cylinder

The CFD solutions were saved every 0.25 s for 500 s resulting in 2000 snapshots. The domain was split into four subdomains, each spanning the entire height of the domain and a quarter of its length. These were discretized with 20 × 20 structured grids. The velocity solutions from the unstructured mesh were linearly interpolated onto the four grids using the finite element basis functions, resulting in a dataset of 8000 samples. For POD, the columns of the snapshots matrix consisted of values of both velocity components, and for the autoencoders, the two velocity components were fed into two separate channels. The training data (which also include the validation data) were formed by randomly selecting 7200 samples from the full dataset. The remaining 800 samples were used as the test dataset (i.e., unseen data).

To test the methods, the solutions are compressed and reconstructed using Eq. (2) for POD, Eq. (6) for the convolutional and adversarial autoencoders, and Eq. (9) for the SVD autoencoder. The error in the reconstruction, Eq. (3), is calculated using the test dataset. Figure 3 shows the effect of the four compression methods (POD, CAE, AAE, and SVD-AE) on a snapshot taken at the 200th time level compared against the original snapshot. The pointwise errors in the velocity magnitude are shown on the right. It can be seen that all four methods (including POD) perform well in their reconstruction of flow past a cylinder. The pointwise errors indicate that, for this snapshot, the convolutional autoencoder gives the best results, followed by the adversarial autoencoder, the SVD-autoencoder, and finally POD. Table I shows the mean of the square reconstruction errors calculated over the test dataset for the flow past a cylinder test case. As seen for the single snapshot in Fig. 3, every compression method that involves an autoencoder outperforms POD.

#### 2. Multiphase flow in a pipe

The domain is split into ten subdomains [each spanning one tenth of the length (*x*) of the domain, but spanning the entire width (*y*) and height (*z*)], which are discretized with 60 × 20 × 20 structured grids. The velocity and volume fraction solutions from the unstructured mesh are interpolated onto these grids over 800 time levels each corresponding to 0.01 s. As before, the finite element basis functions are used to perform the interpolation. The dataset for this test case therefore has a total of 8000 samples (ten subdomains and 800 time levels).

The four compression methods (POD, CAE, AAE, and SVD-AE) are applied to the multiphase flow dataset. For POD, one column of the snapshot matrix consists of nodal values of the three velocity components (each scaled between −1 and 1) and the volume fractions (within the interval [0,1]). For the autoencoders, four channels are used, and scaling is applied to the fields as usual. Initially, ten subdomains were used; however, as the autoencoders were found to have a relatively high error, a further ten subdomains were created that were randomly located within the domain, making a total of 20 subdomains (and 16 000 samples in the dataset). Having more subdomains provided more training data, which probably led to the observed improvement in the results. The autoencoders were trained with 90% of the data chosen at random from the dataset. For details of the hyperparameter optimsation and the networks used, see Tables IV–,VI in the Appendix.

Figure 4 shows how the autoencoders performed in reconstructing the pipe flow dataset. It is not surprising that they seem to perform less well than for the flow past a cylinder case, given the fact that the compression ratio was 80 for flow past a cylinder, whereas, for pipe flow, it was 9600. (For the former a 20 × 20 grid with two fields was compressed to ten variables, whereas for the latter this a 60 × 20 × 20 grid with four fields was compressed to ten variables.) Even at this compression ratio, all dimensionality reduction methods seemed able to reconstruct the slug in Fig. 4 to some degree, with the convolutional AE doing this particularly well. For easier visualization, Fig. 4 shows just part of the domain, which includes a slug and also two boundaries between subdomains. The boundary at 2 m can be identified by a slight kink that can be observed particularly well in the reconstructions of the AAE and the SVD-AE. This kink appears to the left of the slug, and highlights that for some models these boundaries induced additional inaccuracies. This issue could be addressed in future research by allowing the compressive methods to see the solutions of the neighboring subdomains during compression, so that they can explicitly take this boundary into account.

Table II shows the reconstruction error over the test data for the dimensionality reduction methods. Here, Eq. (3) was used, where vectors $uk$ and $(urecon)k$ consist of the scaled velocities and volume fractions. Once again, the convolutional autoencoder has the lowest errors.

### C. Prediction for multiphase flow in a pipe

As the convolutional autoencoder performed better than the other networks for dimensionality reduction, we go on to combine this with a predictive adversarial network within a domain decomposition framework to form a reduced-order model (AI-DDNIROM). A schematic diagram of how the networks are combined can be seen in Fig. 2.

#### 1. Training and predicting with the original domain

For the prediction, the HFM produced 1400 solutions over 14 s of real time. The training and validation data were taken from time levels 1 to 799, and the test data from time levels 800 to 1400. Hyperparameter optimization was performed, and the results of this can be found in Tables VII and VIII of the Appendix. As part of this process, it was found that the best time step for the NIROM was 0.06 s, i.e., 6 times as large as the time interval between the HFM solutions. The MSE achieved on the validation data was $15.6\xd710\u22124$ and on the test data was $103\xd710\u22124$. Figure 5 compares the predictions of volume fraction with those of the HFM and shows the pointwise error for two snapshots in the test data (unseen by the model). The agreement between the predictive adversarial network and the HFM is very good.

#### 2. Extending the domain and associated boundary conditions

Having trained an AI-DDNIROM in Sec. III C I with snapshots from the 10 m long pipe and made predictions for that pipe, in this section. we use the method described in Sec. II D to predict the flow evolution and volume fractions along a pipe of length 98 m based on training data from the 10 m pipe. The extended pipe is split into 98 subdomamins, for which the initial conditions come from the simulation of the 10 m pipe taken at 7.2 s (time level 720). This is in order to start simulating from a state that is well developed. The first subdomain of the 98 m pipe takes initial conditions from the third subdomain of the 10 m pipe; the second to seventh subdomains of the 98 m pipe take the values from the fourth to the ninth subdomains of the 10 m pipe. This is repeated 15 more times, and the final subdomain of the 98 m pipe takes the tenth and final subdomain of the 10 m pipe, see Fig. 6. The first, second, and tenth subdomains of the 10 m pipe were not used, to avoid introducing any spurious effects from the boundaries.

Velocity and volume fractions are specified throughout time in the first and last (98th) subdomains, which act in a manner similar to boundary conditions. There is no high-fidelity model for the 98 m pipe from which to take boundary conditions, and, as the time over which predictions are made exceeds the time over which snapshots were collected from the high-fidelity model of the 10 m pipe, boundary conditions must be generated somehow. Three methods of producing boundary conditions are reported (as described in Sec. II D):

- (i)Cycling through slug formation: a slug is found in the shorter pipe, and the velocity and volume fraction fields associated with the advection of this slug through the subdomain are repeated as required. The particular subdomain of the shorter pipe was the third (between 2 and 3 m), between time levels 750 and 804. So, the boundary condition for the left-most end of the extended pipe can be written as(21)$\alpha 1ext(tk)=\alpha 3(t\u0303k)\u2003\u2200tk\u2009\u2a7e\u20090,$(22)$u1ext(tk)=u3(t\u0303k)\u2003\u2200tk\u2009\u2a7e\u20090,$where $tk=k\Delta t$ for time level
*k*, ($k=0,1,\u2026$,) and a time step of $\Delta t$, and(23)$t\u0303k=(tk\Delta t(mod54)+750)\Delta t.$where $a(modn)$ gives the non-negative remainder when

*n*has been subtracted from*a*as many times as possible. For this example, the time step of the reduced-order model is 0.06 s. The slug appears in this subdomain shortly after the selected time window as a relatively thin instability, of the order of magnitude of 10 cm in length, and develops in width as it advects through the domain. - (ii)Perturbed instability: at the 798th time level an instability occurs in the third subdomain of the shorter pipe. The volume fraction field associated with this is perturbed spatially by Gaussian noise, the velocity field is left unperturbed, and both are used as boundary conditions in the first subdomain of the extended pipe. In the following, $\alpha 1ext$ is the volume fraction in the first subdomain of the extended pipe, $\alpha 3$ is the volume fraction in the third subdomain of the shorter pipe:(24)$\alpha 1ext(tk)=\alpha 3(t\u0303k)+r\u2003\u2200tk\u2009\u2a7e\u20090,$(25)$u1ext(tk)=u3(t\u0303k)\u2003\u2200tk\u2009\u2a7e\u20090,$
where $t\u0303k=$ 7.98 s and

*r*is a random spatial perturbation. - (iii)Original boundaries repeated: velocity and volume fraction solutions from the third subdomain of the shorter pipe are used as the boundary conditions for the first subdomain in the extended pipe and repeated for as long as required. The solution fields from time 1 s to 8 s are used, as this corresponds to times where the air had passed through the entire length of the shorter pipe. Therefore, for times in $[0,7)$ seconds of the extended pipe, times in $[1,8)$ seconds of the shorter pipe are used; for times in [7,14) seconds of the extended pipe, times in [1,8) seconds of the shorter pipe are used; etc. So, the boundary condition for the left-most end of the extended pipe can be written as(26)$\alpha 1ext(tk)=\alpha 3(t\u0303k)\u2003\u2200tk\u2009\u2a7e\u20090,$(27)$u1ext(tk)=u3(t\u0303k)\u2003\u2200tk\u2009\u2a7e\u20090,$where $tk=k\Delta t$ for time level
*k*, ($k=0,1,\u2026$,) and a time step of $\Delta t$, and(28)$t\u0303k=(tk\Delta t(mod700)+100)\Delta t.$

In all cases, the boundary condition for the final subdomain is based on a snapshot from subdomain 2 at 7.5 s in the 10 m pipe when the flow was almost steady with the lower half of the pipe occupied by water and the upper half occupied by air.

Various statistics are presented in this section in an aim to assess whether the AI-DDNIROM approach produces realistic results. If the expected advantage of the adversarial training strategy to produce a model that does not extrapolate beyond the seen training data holds true, then the predictive model could be expected to not diverge significantly from the original simulation. Figures 7(a)–7(c) show how the liquid volume fraction field varies over time in the original simulation and for the reduced models using two of the tested boundary conditions (cycling through slug formation and perturbed instability). The results obtained when repeating the original boundary conditions were similar to the original simulation and are not shown here. The time interval for the two reduced-order models corresponds to the instabilities having passed through two thirds of the pipe. Time series data were collected at values of $x=6.5\u2009m$ (for the original simulation) and $x=64.5\u2009m$ (for the two reduced models), and at a height of 0.0039 m (a tenth of the pipe radius) above the centerline of the pipe. To analyze the frequency spectra, a discrete Fourier transform was then applied to the data. Figures 7(d)–7(f) show that the slug characteristic frequency spectra for the predictions are similar to that of the original simulation. In particular, the main peak has a similar value in all three simulations (original simulation: 0.76 Hz; reduced model which cycles through slug formation: 0.7 Hz; reduced model with the perturbed instability 0.88 Hz). This suggests that simulations from the AI-DDNIROMs based on either of these boundary conditions are able to behave in a realistic way. In fact, the frequency of the main peak could be interpreted as the pseudo-slug frequency. Technically slugs are only defined as such when they span the full vertical extent of the pipe. On the other hand, pseudo-slugs^{90} or proto-slugs^{91} are precursors to slugs, which do not necessarily reach the full height of the pipe.

Figure 8 follows a pseudo-slug for five time levels as viewed through the volume fraction fields for the original simulation and the reduced-order models. It shows first that the instabilities presented in Figs. 8(d) and 8(e) were similar to an instability that also occurred within the original simulation, presented in Fig. 8(c). Furthermore, by observing that the instabilities traveled similar distances between time levels, it can be deduced that they traveled at a similar velocity within the shown timespan as well. While this only presents the dynamics for a single instability at a couple of points in time, the similarity of these situations might reveal that the predictive adversarial model was producing a situation very similar to one it had seen before, which is what this model was hypothesized to do due to its use of the adversarial training strategy.

Figures 9(a)–9(c) show the volume fractions averaged over the full 98 m domain for a short time period. Note that the start of this time period was chosen so that the influence of the boundaries had already propagated throughout the domain. Figures 9(d)–9(f) display the volume fractions averaged over the time period included in the previous three subfigures [Figs. 9(a)–9(c)], and the width and length of the domain. These latter three plots thus show how the volume fractions change with respect to height. It is clear from these plots that most of the water collects at the bottom of the pipe. If we assume that the situation in which the original boundaries were repeated in their entirety produced results similar to the original simulation, the fact that the plotted dynamics are similar to the two simulations with artificially generated boundaries suggests strongly that the model was making realistic predictions here for each of the boundaries.

Figure 10 displays the volume fractions predicted by the AI-DDNIROMs throughout the pipe along the centerline (*xz* plane) for 12 s in time and averaged over the height of the pipe. The red bands that stretch diagonally along the plotted domain represent slugs propagating downstream through the domain as time progresses. The slopes of these lines represent the corresponding velocities. Black lines have been drawn on these plots to indicate these liquid slug velocities and also velocities of secondary waves (light blue). The velocity magnitudes are given in Table III. From these plots and the table, one can see that both the slug velocities and velocities of the secondary waves produced by different boundary conditions are very similar. Note that the slight variations may have been caused by the interactions of slugs being within each others vicinity. A pattern that is generally observed in each of the three graphs in Fig. 10 is that two slugs which are close to one another tend to slowly approach one another. The slug which is ahead seems to slow down and disappear as the approaching slug catches up. These graphs also clearly display the influence of the boundary at the inlet on the simulation. In fact, the first couple of meters is where the simulations differ the most. However, it seems that after those first couple of meters the simulations all restore to a very similar pattern. In a similar experimental setup to the computational domain modeled here, what has been observed is that the slug frequencies are largest near the inlet, but, after about 5 m, settle to a value independent of the distance from the inlet.^{92} Correlations for slug frequency are often sensitive to the superficial liquid velocity, which, in turn depends on the mass of liquid in the pipe,^{93} so the fact that the slug frequencies do not appear to change significantly after about 10 m (see corresponding slug lengths in Fig. 10), seems to suggest that the mass of liquid is conserved along the pipe.

. | Slug velocity . | Secondary wave velocity . |
---|---|---|

Cycling through slug formation | 3.5 m s^{−1} | 1.5 m s^{−1} |

Perturbed instability | 3.4 m s^{−1} | 1. m s^{−1} |

Original boundaries repeated | 3.3 m s^{−1} | 1.5 m s^{−1} |

. | Slug velocity . | Secondary wave velocity . |
---|---|---|

Cycling through slug formation | 3.5 m s^{−1} | 1.5 m s^{−1} |

Perturbed instability | 3.4 m s^{−1} | 1. m s^{−1} |

Original boundaries repeated | 3.3 m s^{−1} | 1.5 m s^{−1} |

#### 3. Computational times

The AI-DDNIROM shows a significant computational speed up over the high-fidelity model as expected. The high-fidelity model of the 10 m pipe took approximately two weeks to complete (run on processor type Intel^{®} Xeon^{®} E5–2640, 2.4 GHz), whereas the AA-DDNIROM prediction for the 98 m pipe took approximately 20 minutes to generate (run on GPUs within Google's Colab platform^{94,95})

## IV. CONCLUSIONS AND FURTHER WORK

We present an AI-based non-intrusive reduced-order model combined with domain decomposition (AI-DDNIROM), which is capable of making predictions for significantly larger domains than the domain used in training the model. For dimensionality reduction we use a convolutional autoencoder and for prediction we use a predictive adversarial network. During training, the predictive adversarial network learns the underlying probability distribution (of the latent variables) associated with the fluids' behavior. The main findings of this study are listed below.

- (i)
For dimensionality reduction, a number of autoencoders are compared with proper orthogonal decomposition, and the convolutional autoencoder is seen to perform the best for both test cases (2D flow past a cylinder and 3D multiphase flow in a pipe).

- (ii)
When training neural networks, it has been observed that computational physics applications typically have access to less training data than image-based applications,

^{12,23}which can lead to poor generalization. To combat this, for the dimensionality reduction of multiphase flow in a pipe, we use “overlapping” snapshots, that is, in addition to ten subdomains being equally spaced along the pipe, ten supplementary subdomains are located at random within the pipe. This doubles the amount of training data and results in improved performance. - (iii)
For prediction, we use an predictive adversarial network based on the adversarial autoencoder

^{53}but modified to predict in time. This model performs well, gives realistic results, and, unlike feed forward or recurrent networks without such an adversarial layer, does not diverge for the multiphase test case shown here. - (iv)
Finally, we make predictions for a 98 m pipe (the “extended pipe”) with the AI-DDNIROM that was trained on results from a 10 m pipe. Statistics of results from the extended pipe are similar to those of the original pipe, so we conclude that the predictive adversarial network has made realistic predictions for this extended pipe.

A number of improvements could be made to the approach presented here. A physics-informed term could be included in the loss function of either the convolutional autoencoder or the predictive adversarial network. This would ensure that conservation of mass and momentum would be more closely satisfied by the predictions of the neural networks. Second, although the initial conditions have little effect on the predictions, the boundary conditions do have a significant effect. Rather than the heuristic approach adopted here, a generative adversarial model (GAN) could be used to predict boundary conditions for the inlet and outlet subdomains. The GAN could be trained to predict the reduced variables at several time levels, then latent variables consistent with all but one of the time levels (the future time level) can be found by an optimization approach.^{67} From these latent variables, the boundary condition for the future time level can be obtained. Finally, a hierarchy of reduced-order models could be used in order to make the approach faster. The lowest-order model could represent the simplest physical features of the flow, and the higher-order models could represent more complicated flow features. To decide whether the model being used in a particular subdomain was sufficient, the discriminator of the predictive adversarial network could be used.

## ACKNOWLEDGMENTS

The authors would like to acknowledge the following EPSRC grants: MUFFINS, MUltiphase Flow-induced Fluid-flexible structure InteractioN in Subsea applications (Nos. EP/P033180/1 and EP/P033148/1); RELIANT, Risk EvaLuatIon fAst iNtelligent Tool for COVID19 (No. EP/V036777/1); the PREMIERE programme grant (No. EP/T000414/1); MAGIC, Managing Air for Green Inner Cities (No. EP/N010221/1); and INHALE, Health assessment across biological length scales (No. EP/T003189/1). We would also like to acknowledge the Applied Computational Science and Engineering MSc course at Imperial, during which, Zef Wolffs obtained the results shown in this paper. Finally, we would like to thank the reviewers for their comments and suggestions, which have improved the paper.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

C.E.H.: conceptualization, methodology, software, writing—original draft, writing—review and editing, supervision; Z.W.: methodology, software, writing —original draft, writing—review and editing; J.A.T.: software, writing—review and editing; L.K.: software, writing—review and editing; P.S.: software, writing—review and editing; A.N. conceptualization, writing—review and editing, supervision; I.M.N.: conceptualization, writing—review and editing; O.K.M.: conceptualization, writing—review and editing, supervision, funding acquisition; N.S.: conceptualization, writing—original draft, writing—review and editing, supervision, funding acquisition; C.C.P.: conceptualization, methodology, software, writing—original draft, writing—review and editing, supervision, funding acquisition.

## DATA AVAILABILITY

The data that support the findings of this study are openly available in Github at https://github.com/acse-zrw20/DD-GAN-AE, Ref. 96.

### APPENDIX: HYPERPARAMETER OPTIMIZATION

Extensive hyperparameter optimzation was carried out for the artificial neural networks used in this investigation. This was done on the Weights and Biases platform which allows for efficient searching of high-dimensional parameter space, using methods such as random searches and Bayesian searches. For example, to perform a grid search of the predictive adversarial network for one architecture would involve searching 18 dimensional parameter space, and, with the combinations given in Table IV, would amount to over 2 billion ($2\xd7109$) model evaluations (for one architecture). Instead of using a grid search, we perform an initial random search of parameter space, followed by Bayesian optimization. For the predictive adversarial network, this resulted in 1530 model evaluations (for all architectures). The full report for this network is available on Weights & Biases.

All networks . | |
---|---|

Activation functions | tanh, sigmoid, relu, elu |

Final activation Function | tanh, sigmoid, linear |

Architecture^{a} | Number of layers: 6, $\cdots $, 20 |

Number of channels: 2, $\cdots $, 128 | |

Dense layer sizes (non-latent): 32, $\cdots $, 2000 | |

Kernel sizes: 3, 5 | |

Layer types: {1D, 2D, 3D}-Conv., {1D, 2D, 3D}-MaxPool, {1D, 2D, 3D}-UpSample, Dense | |

Batch size | 32, 64, 128 |

Optimizer | Adam, Nadam, SGD |

β_{1} | 0.8, 0.9, 0.98 |

β_{2} | 0.9, 0.999, 1 |

Batch normalization | True, false |

Dropout | 0.3, 0.55, 0.8 |

Epochs | 100, 200, 500, 1000, 2000 |

Interval | 1, 2, 4, 5, 6, 10 |

Learning rate | 0.000 05, 0.0005, 0.005 |

Adversarial networks only | |

Discrim architecture^{a} | Number of layers: 3 |

Dense layer sizes (non-latent): 100, 500, 1000 | |

Layer types: dense | |

n discrim | 1, 2, 5 |

n gradient | 0, 3, 8, 15, 30 (0 means that no steps of gradient ascent were taken) |

std noise | 0, 0.000 01, 0.001, 0.01, 0.05, 0.1 |

Regularization | 0, 0.000 001, 0.000 01, 0.001 |

Predictive adversarial networks only | |

Latent vars | 30, 50, 100 |

All networks . | |
---|---|

Activation functions | tanh, sigmoid, relu, elu |

Final activation Function | tanh, sigmoid, linear |

Architecture^{a} | Number of layers: 6, $\cdots $, 20 |

Number of channels: 2, $\cdots $, 128 | |

Dense layer sizes (non-latent): 32, $\cdots $, 2000 | |

Kernel sizes: 3, 5 | |

Layer types: {1D, 2D, 3D}-Conv., {1D, 2D, 3D}-MaxPool, {1D, 2D, 3D}-UpSample, Dense | |

Batch size | 32, 64, 128 |

Optimizer | Adam, Nadam, SGD |

β_{1} | 0.8, 0.9, 0.98 |

β_{2} | 0.9, 0.999, 1 |

Batch normalization | True, false |

Dropout | 0.3, 0.55, 0.8 |

Epochs | 100, 200, 500, 1000, 2000 |

Interval | 1, 2, 4, 5, 6, 10 |

Learning rate | 0.000 05, 0.0005, 0.005 |

Adversarial networks only | |

Discrim architecture^{a} | Number of layers: 3 |

Dense layer sizes (non-latent): 100, 500, 1000 | |

Layer types: dense | |

n discrim | 1, 2, 5 |

n gradient | 0, 3, 8, 15, 30 (0 means that no steps of gradient ascent were taken) |

std noise | 0, 0.000 01, 0.001, 0.01, 0.05, 0.1 |

Regularization | 0, 0.000 001, 0.000 01, 0.001 |

Predictive adversarial networks only | |

Latent vars | 30, 50, 100 |

^{a}

Here a global picture of the architectures is presented, for the source code containing all of the used architectures please see the Github repository: https://github.com/acse-zrw20/DD-GAN-AE/tree/main/ddganAE/architectures.

Table IV shows the range of hyperparameters that were investigated during optimization for all the networks (the three autoencoder-based networks used for dimensionality reduction for the two test cases and the predictive adversarial network used in multiphase flow in a pipe). These include the exponential decay rate for the first moment estimates (*β*_{1}); the exponential decay rate for the exponentially weighted infinity norm (*β*_{2}); the interval between snapshots (interval) so that an interval of *n* corresponds to every *n*th snapshot being put in the datasets; the number of discriminator iterations (n discrim); the number of gradient ascent steps (n gradient); the standard deviation of the noise that was randomly added to the input of the discriminator within the adversarial autoencoder.

Table V shows the optimal values found in the hyperparameter optimization for the dimensionality reduction methods based on autoencoders for flow past a cylinder and for multiphase flow in a pipe.

. | Flow past a cylinder . | Multiphase pipe flow . | ||||
---|---|---|---|---|---|---|

. | CAE . | AAE . | SVD-AE . | CAE . | AAE . | SVD-AE . |

Activation functions: | ||||||

Convolutional layers | elu | elu | elu | elu | sigmoid | |

Dense layers | relu | relu | relu | relu | linear | sigmoid |

Output layer | elu | sigmoid | sigmoid | linear | ||

Optimizer: | ||||||

Method | Adam | Nadam | Nadam | Adam | Adam | Nadam |

β_{1} | 0.98 | 0.9 | 0.98 | 0.8 | 0.9 | 0.8 |

β_{2} | 0.9 | 0.999 99 | 0.999 99 | 0.9 | 0.9 | 0.999 99 |

Batch size | 128 | 128 | 64 | 64 | 32 | 64 |

Epochs | 200 | 200 | 200 | 100 | 1000 | 100 |

Batch normalization | False | False | ||||

Train method | Default | Default | ||||

Dropout | 0.55 | 0.55 | ||||

Learning rate | 0.00005 | 0.000 005 | 0.0005 | 0.0005 | 0.000 005 | 0.000 05 |

Regularization | 0 | 0 | 0 | 0 | 0 | 0 |

. | Flow past a cylinder . | Multiphase pipe flow . | ||||
---|---|---|---|---|---|---|

. | CAE . | AAE . | SVD-AE . | CAE . | AAE . | SVD-AE . |

Activation functions: | ||||||

Convolutional layers | elu | elu | elu | elu | sigmoid | |

Dense layers | relu | relu | relu | relu | linear | sigmoid |

Output layer | elu | sigmoid | sigmoid | linear | ||

Optimizer: | ||||||

Method | Adam | Nadam | Nadam | Adam | Adam | Nadam |

β_{1} | 0.98 | 0.9 | 0.98 | 0.8 | 0.9 | 0.8 |

β_{2} | 0.9 | 0.999 99 | 0.999 99 | 0.9 | 0.9 | 0.999 99 |

Batch size | 128 | 128 | 64 | 64 | 32 | 64 |

Epochs | 200 | 200 | 200 | 100 | 1000 | 100 |

Batch normalization | False | False | ||||

Train method | Default | Default | ||||

Dropout | 0.55 | 0.55 | ||||

Learning rate | 0.00005 | 0.000 005 | 0.0005 | 0.0005 | 0.000 005 | 0.000 05 |

Regularization | 0 | 0 | 0 | 0 | 0 | 0 |

Table VI gives the optimal architectures found by hyperparameter optimization the six autoencoder-based networks used in the dimensionality reduction of flow past a cylinder and multiphase flow in a pipe.

. | Flow past a cylinder . | Multiphase flow in a pipe . | ||||
---|---|---|---|---|---|---|

Layers . | CAE . | AAE . | SVD-AE . | CAE . | AAE . | SVD-AE . |

Input | (55, 42, 2) | (55, 42, 2) | 100 | (60, 20, 20, 4) | (60, 20, 20, 4) | 100 |

Conv | (55, 42, 32) | (55, 42, 32) | (60, 20, 20, 32) | (60, 20, 20, 32) | ||

MaxPool | (28, 21, 32) | (28, 21, 32) | (30, 10, 10, 32) | (30, 10, 10, 32) | ||

Conv | (28, 21, 64) | (28, 21, 64) | (30, 10, 10, 64) | (30, 10, 10, 64) | ||

MaxPool | (14, 11, 64) | (14, 11, 64) | (15, 5, 5, 64) | (15, 5, 5, 64) | ||

Conv | (14, 11, 128) | (15, 5, 5, 128) | (15, 5, 5, 128) | |||

MaxPool | (7, 6, 128) | (8, 3, 3, 128) | (8, 3, 3, 128) | |||

Flatten | 5376 | 9856 | 9216 | 9216 | ||

Dense 1 | 2688 | 9856 | 500^{a} | 10 | 4608 | 1500 |

Dense 2 | 10 | 4926 | 500^{a} | 9216 | 10 | 2000 |

Dense 3 | 2688 | 10 | 10 | 4608 | 10 | |

Dense 4 | 5376 | 4926 | 500^{a} | 9218 | 1500 | |

Dense 5 | 9856 | 500^{a} | 2000 | |||

Dense 6 | 9856 | |||||

Reshape | (7, 6, 128) | (14, 11, 64) | (8, 3, 3, 128) | (8, 3, 3, 128) | ||

Conv | (7, 6, 128) | (14, 11, 64) | (8, 3, 3, 128) | (8, 3, 3, 128) | ||

UpSample | (14, 12, 128) | (28, 22, 64) | (16, 6, 6, 128) | (16, 6, 6, 128) | ||

Conv | (14, 12, 64) | (28, 22, 32) | (16, 6, 6, 64) | (16, 6, 6, 64) | ||

Upsample | (28, 24, 64) | (56, 44, 32) | (32, 12, 12, 64) | (32, 12, 12, 64) | ||

Conv | (28, 24, 32) | (56, 44, 2) | (30, 10, 10, 32)^{b} | (30, 10, 10, 32)^{b} | ||

UpSample | (56, 48, 32) | (60, 20, 20, 32) | (60, 20, 20, 32) | |||

Conv | (56, 48, 2) | (60, 20, 20, 4) | (60, 20, 20, 4) | |||

Crop | (55, 42, 2) | |||||

Output | 100 | 100 | ||||

Trainable parameters | 29 300 300 | 291 587 010 | 612 110 | 1196 238 | 86 001 860 | 6 392 110 |

. | Flow past a cylinder . | Multiphase flow in a pipe . | ||||
---|---|---|---|---|---|---|

Layers . | CAE . | AAE . | SVD-AE . | CAE . | AAE . | SVD-AE . |

Input | (55, 42, 2) | (55, 42, 2) | 100 | (60, 20, 20, 4) | (60, 20, 20, 4) | 100 |

Conv | (55, 42, 32) | (55, 42, 32) | (60, 20, 20, 32) | (60, 20, 20, 32) | ||

MaxPool | (28, 21, 32) | (28, 21, 32) | (30, 10, 10, 32) | (30, 10, 10, 32) | ||

Conv | (28, 21, 64) | (28, 21, 64) | (30, 10, 10, 64) | (30, 10, 10, 64) | ||

MaxPool | (14, 11, 64) | (14, 11, 64) | (15, 5, 5, 64) | (15, 5, 5, 64) | ||

Conv | (14, 11, 128) | (15, 5, 5, 128) | (15, 5, 5, 128) | |||

MaxPool | (7, 6, 128) | (8, 3, 3, 128) | (8, 3, 3, 128) | |||

Flatten | 5376 | 9856 | 9216 | 9216 | ||

Dense 1 | 2688 | 9856 | 500^{a} | 10 | 4608 | 1500 |

Dense 2 | 10 | 4926 | 500^{a} | 9216 | 10 | 2000 |

Dense 3 | 2688 | 10 | 10 | 4608 | 10 | |

Dense 4 | 5376 | 4926 | 500^{a} | 9218 | 1500 | |

Dense 5 | 9856 | 500^{a} | 2000 | |||

Dense 6 | 9856 | |||||

Reshape | (7, 6, 128) | (14, 11, 64) | (8, 3, 3, 128) | (8, 3, 3, 128) | ||

Conv | (7, 6, 128) | (14, 11, 64) | (8, 3, 3, 128) | (8, 3, 3, 128) | ||

UpSample | (14, 12, 128) | (28, 22, 64) | (16, 6, 6, 128) | (16, 6, 6, 128) | ||

Conv | (14, 12, 64) | (28, 22, 32) | (16, 6, 6, 64) | (16, 6, 6, 64) | ||

Upsample | (28, 24, 64) | (56, 44, 32) | (32, 12, 12, 64) | (32, 12, 12, 64) | ||

Conv | (28, 24, 32) | (56, 44, 2) | (30, 10, 10, 32)^{b} | (30, 10, 10, 32)^{b} | ||

UpSample | (56, 48, 32) | (60, 20, 20, 32) | (60, 20, 20, 32) | |||

Conv | (56, 48, 2) | (60, 20, 20, 4) | (60, 20, 20, 4) | |||

Crop | (55, 42, 2) | |||||

Output | 100 | 100 | ||||

Trainable parameters | 29 300 300 | 291 587 010 | 612 110 | 1196 238 | 86 001 860 | 6 392 110 |

^{a}

Denotes a layer which is followed by a dropout layer during training.

^{b}

Denotes convolutional layers which have no padding. In all other cases padding is set so that the output has the same dimensions as the input array, although the number of channels may vary.

Table VII shows the optimal values found in the hyperparameter optimization for the predictive adversarial network used for the non-intrusive reduced-order model of multiphase flow in a pipe, and Table VIII gives the optimal architecture.

Activation functions: | Dropout | 0.3 | |

Convolutional layers | relu | Interval | 6 |

Dense layers | relu | Learning rate | 0.000 05 |

Final layer | tanh | Latent vars | 100 |

Optimizer: | n discrim | 1 | |

Method | Nadam | n gradient | 15 |

β_{1} | 0.98 | std noise | 0.01 |

β_{2} | 0.9 | Regularization | 0.001 |

Batch size | 32 | Batch normalization | True |

Epochs | 2000 | Training method | Weighted loss |

Activation functions: | Dropout | 0.3 | |

Convolutional layers | relu | Interval | 6 |

Dense layers | relu | Learning rate | 0.000 05 |

Final layer | tanh | Latent vars | 100 |

Optimizer: | n discrim | 1 | |

Method | Nadam | n gradient | 15 |

β_{1} | 0.98 | std noise | 0.01 |

β_{2} | 0.9 | Regularization | 0.001 |

Batch size | 32 | Batch normalization | True |

Epochs | 2000 | Training method | Weighted loss |

Layers . | Predictive adversarial network . | Discriminator . |
---|---|---|

Input | 30 | 100 |

Dense 1 | 500 | 100 |

Dense 2 | 500 | 500 |

Dense 3 | 100 | |

Dense 4 | 500 | |

Dense 5 | 500 | |

Output | 10 | 1 |

Trainable parameters | 622 110 | 61 101 |

Layers . | Predictive adversarial network . | Discriminator . |
---|---|---|

Input | 30 | 100 |

Dense 1 | 500 | 100 |

Dense 2 | 500 | 500 |

Dense 3 | 100 | |

Dense 4 | 500 | |

Dense 5 | 500 | |

Output | 10 | 1 |

Trainable parameters | 622 110 | 61 101 |