An AI-based Domain-Decomposition Non-Intrusive Reduced-Order Model for Extended Domains applied to Multiphase Flow in Pipes

The modelling 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 kilometres versus a pipe diameter of just a few inches. 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 using domain decomposition; dimensionality reduction; training a neural network to make predictions for a single subdomain; and by using an iteration-by-subdomain technique to converge the solution over the whole domain. To find the low-dimensional space, we explore several types of autoencoder networks, known for their ability to compress information accurately and compactly. The performance of the autoencoders is assessed on 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 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. Statistics of the flows obtained from the CFD simulations are compared to those of the AI-DDNIROM predictions to demonstrate the success of our approach.

efficient than a linear subspace for approximating the HFM.
Convolutional networks are particularly good at analysing 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][34][35]; the Sod shock test and transient wake of a ship [36]; air pollution in an urban environment [37][38][39]; parametrised 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 parametrised 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][42][43][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 high Reynolds number [49] (described as hybrid ROM); the generalised 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 parametrised 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 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][55][56][57][58][59][60], Gaussian Process Regression [11,45,[61][62][63] and Long-Short Term Memory networks [31,32,34,35,38,51,64]. When using these types of neural network 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 References [56,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 modelling. Two exceptions are Reference [66] which combines a VAE and GAN to model flow past a cylinder and the collapse of a water dam; and Reference [67] which uses a GAN to predict the reduced variables of an epidemiological model which modelled the spread of a virus through a small, idealised 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 hydromechanical 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 nonintrusive ROMs [61], which was later extended to partition the domain by minimising 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 behaviour 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 Figure 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 Reference [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 modelling 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-tophase 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][79][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, the phase interfacial interactions), depending on the flow pattern, flow direction and pipe physical properties (inclination, diameter, length).
Significant progress has been made in 3D modelling [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 advectiondominated 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. 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; (iii) the an 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 2 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 3 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 optimisation process and the network architectures are given in the appendix.

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 discretised governing equations (the high-resolution or high-fidelity model); (ii) reducing the dimensionality of the discretised 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 (Section 2.2) and prediction (Section 2.3) are described. The final section (Section 2.4) outlines an approach for making predictions for a larger domain having used a smaller domain to generate the training data.

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.

Proper Orthogonal Decomposition
Proper Orthogonal Decomposition is a commonly used technique for dimensionality reduction when constructing reduced-order models. POD requires the minimisation 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 minimise 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 S can be decomposed as where the matrix U contains the left singular vectors, V 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 U can be discarded. The POD basis functions to be retained are stored in a matrix denoted by R. The POD coefficients of a snapshot can be found by pre-multiplying the snapshot by R T , and the reconstruction of a snapshot can be found by pre-multiplying the POD coefficients of the snapshot by R: where u k is the kth snapshot and (u recon ) k is its reconstruction. Hence the reconstruction error over a set of N snapshots {u 1 , u 2 , . . . , u N } 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 snapshots matrix: where u i and v i represent the x and y components of velocity, respectively, at the ith 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 w k i and α k i represent the z component of velocity and the volume fraction, respectively, at the ith node of the kth snapshot. In this case, the velocity components are scaled to be in the range [−1, 1] so that their magnitudes are similar to those of the volume fractions.

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 u k represents the kth sample in the dataset of N samples and (u recon ) k represents the corresponding output of the autoencoder, which can be written as then the mean square error can be expressed as in Equation (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 (P prior ). 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 minimised (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 Equation (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 ∼ P prior is a sample from the desired distribution and u ∼ P data is a sample input grid. There are strong similarities between the adversarial autoencoder, GANs and Variational Autoencoders (VAEs). All three types of network set out to obtain better generalisation than non-adversarial networks by attempting to obtain a smooth latent space with no gaps. Results in Reference [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].

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 f enc is the encoder, R represents the POD basis functions, u k is the kth snapshot and z k 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 f dec is the decoder, f ae is the autoencoder and (u recon ) k is the reconstruction of the kth 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.

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 timestep. The adversarial autoencoder is re-purposed for this task in an attempt to capitalise on the fact that this network should produce realistic results (providing that the training data is representative of the behaviour that will be modelled). 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 neighbouring 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 minimising the error between its output and the predicted variables at the future time level, as well as incorporating the adversarial training strategy described in Section 2.2.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 next section). Given the aspect ratio of the pipe, we split the domain into subdomains of equal length in the axial direction, see where z k i represents the reduced variables in subdomain i at the future time level k; z k−1 i represents the same but at the preceding time level; and z k i−1 and z k i+1 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 (deceasing i). During the iteration process, z k i of Equation (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 z k i−1 and z k i+1 in Equation (10)). The solution at the previous, converged time level could be used (z k i±1 = z k−1 i±1 ), 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 if k > 1 then 8: for subdomain i = 2, 3, . . . , N sub − 1 do 9: end for !! sweep over subdomains 17: for sweep iteration j = 1, 2, . . . , N sweep do 18: !! calculate the latent variables of subdomain i at time level k 20: end for 22: end for 23: end for

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 Section 2.3).
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) was 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.

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 test cases are now described.

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, τ contains the viscous terms associated with an isotropic Newtonian fluid, p represents the non-hydrostatic pressure, t is time and the gradient operator ∇ is defined as When solving these equations, a linear triangular element is used with a discontinuous Galerkin discretisation for the velocities and a continuous Galerkin representation of the pressure (the P1DG-P1 element). To discretise in time, Crank-Nicolson is used. 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 discretised and solved, the reader is referred to [88].
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 σ is the force representing surface tension. The bulk density and bulk viscosity are defined as respectively, where ρ water and ρ air are the densities of water and air respectively, and µ water and µ air are the dynamic viscosities of water and air respectively. For more details of how the governing equations are discretised and solved, see Reference [86], including information on the unstructured adaptive meshing

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 A.4, A.5 and A.6. 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 Figure 2 shows a schematic diagram of how the networks used for dimensionality reduction are trained for the flow past a cylinder test case.

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 discretised with 20 by 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 includes the validation data) was 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 Equation (2) for POD, Equation (6) for the convolutional and adversarial autoencoders, and Equation (9) for the SVD autoencoder. The error in the reconstruction, Equation (3), is calculated using the test dataset. Figure 3 shows the effect of  Table 1 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 Figure 3, every compression method that involves an autoencoder outperforms POD.

Multiphase flow in a pipe
The domain is split into 10 subdomains (each spanning one tenth of the length (x) of the domain, but spanning the entire width (y) and height (z)), which are discretised with 60 by 20 by 20 structured grids.
The velocity and volume fraction solutions from the unstructured mesh are interpolated onto these grids over channels are used and scaling is applied to the fields as usual. Initially 10 subdomains were used, however, as the autoencoders were found to have a relatively high error, a further 10 subdomains were created that were randomly located within the domain, making a total of 20 subdomains (and 16000 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 A.4, A.5 and A.6 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 by 20 grid with two fields was compressed to 10 variables, whereas for the latter this a 60 by 20 by 20 grid with four fields was compressed to 10 variables.) Even at this compression ratio, all dimensionality reduction methods seemed able to reconstruct the slug in Figure 4 to some degree, with the convolutional AE doing this particularly well. For easier visualisation, Figure 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 neighbouring subdomains during compression, so that they can explicitly take this boundary into account.  Once again, the convolutional autoencoder has the lowest errors.

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  and on the test data was 103 × 10 −4 . 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.

Extending the domain and associated boundary conditions
Having trained an AI-DDNIROM in the previous section with snapshots from the 10 m long pipe and made predictions for that pipe, in this section we use the method described in Section 2.4 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 Figure 6. The first, second and tenth subdomains of the 10 m pipe were not used, to avoid introducing any spurious effects from the boundaries. to 804. So, the boundary condition for the left-most end of the extended pipe can be written as where t k = k∆t for time level k, (k = 0, 1, . . .) and a time step of ∆t, and t k = t k ∆t (mod 54) + 750 ∆t .
where a (mod n) 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, α ext 1 is the volume fraction in the first subdomain of the extended pipe, α 3 is the volume fraction in the third subdomain of the shorter pipe wheret k = 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 where t k = k∆t for time level k, (k = 0, 1, . . .) and a time step of ∆t, and t k = t k ∆t (mod 700) + 100 ∆t .
In all cases, the boundary condition for the final subdomain is based on a snapshot from subdomain 2 at 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 [89] or proto-slugs [90] are precursors to slugs, which do not necessarily reach the full height of the pipe. and 8e were similar to an instability that also occurred within the original simulation, presented in Figure 8c.
Furthermore, by observing that the instabilities travelled similar distances between time levels, it can be deduced that they travelled 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.        Table 3. 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 Figure 10 is that two slugs which are close to one another tend to slowly approach one another. These graphs also clearly display the influence of the boundary at the inlet on the simulation. In fact, the first couple of metres is where the simulations differ the most. However, it seems that after those first couple of metres the simulations all restore to a very similar pattern. Figure 11 shows the volume fractions at a few steps in time to give an impression of these fields throughout the full width of the domain, showing a different perspective of the information seen in Figure 10.  [91,92]).

Conclusions and Further Work
We present an AI-based non-intrusive reduced-order model combined with a domain decomposition (AI-DDNIROM), which is capable of making predictions for significantly larger domains than the domain used in training the model. 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 generalisation. To combat this, for the dimensionality reduction of multiphase flow in a pipe, we use 'overlapping' snapshots, that is, in addition to 10 subdomains being equally spaced along the pipe, 10 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. Secondly, 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 optimisation 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.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request. Some codes and information about the various neural networks used in this paper can be found in the following Github repository: https://github.com/acse-zrw20/DD-GAN-AE.

Appendix A. Hyperparameter optimisation
Extensive hyperparameter optimsation was carried out for the artificial neural networks used in this investigation. This was done on the Weights & 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 A.4, would amount to over 2 billion (2 × 10 9 ) model evaluations (for one architecture). Instead of using a grid search, we perform an initial random search of parameter space, followed by Bayesian optimisation. 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. Table A.4 shows the range of hyperparameters that were investigated during optimisation 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 nth 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 A.5 shows the optimal values found in the hyperparameter optimisation for the dimensionality reduction methods based on autoencoders for flow past a cylinder and for multiphase flow in a pipe. Table A. 6 gives the optimal architectures found by hyperparameter optimisation the six autoencoder-based networks used in the dimensionality reduction of flow past a cylinder and multiphase flow in a pipe. Table A.7 shows the optimal values found in the hyperparameter optimisation for the predictive adversarial network used for the non-intrusive reduced-order model of multiphase flow in a pipe, and      (Conv), maxpooling (MaxPool) or upsampling (UpSample). Flatten layers take an n-dimensional array as an input and return a 1D array as an output. A reshape layer converts a 1D input to have the indicated output dimensions. † denotes a layer which is followed by a dropout layer during training. ‡ 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.