A domain decomposition non-intrusive reduced order model for turbulent flows

Computers and Fluids 182 (2019) 15–27

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

A domain decomposition non-intrusive reduced order model for turbulent flows D. Xiao a,b,c,∗, C.E. Heaney a, F. Fang a,b, L. Mottet a,d, R. Hu a, D.A. Bistrian e, E. Aristodemou a,f, I.M. Navon g, C.C. Pain a,b a

Applied Modelling and Computation Group, Department of Earth Science and Engineering, Imperial College London, Prince Consort Road, London, SW7 2BP, UK Data Assimilation Lab, Data Science Institute, Imperial College London, Prince Consort Road, London, SW7 2BP, UK c ZCCE College of Engineering, Swansea University, Bay Campus, Fabian Way, Swansea, SA1 8EN, UK d Department of Architecture, University of Cambridge, 1-5 Scroope Terrace, Trumpington Street, Cambridge, CB2 IPX, UK e Politehnica University of Timisoara, Hunedoara, 331128, Romania f School of Engineering, London South Bank University, London, UK g Department of Scientific Computing, Florida State University, Tallahassee, FL, 32306-4120, USA b

a r t i c l e

i n f o

Article history: Received 29 August 2018 Revised 25 January 2019 Accepted 14 February 2019 Available online 15 February 2019 Keywords: Non-intrusive reduced order modelling Domain decomposition Machine learning Gaussian process regression Urban flows Turbulent flows Finite element method

a b s t r a c t In this paper, a new Domain Decomposition Non-Intrusive Reduced Order Model (DDNIROM) is developed for turbulent flows. The method works by partitioning the computational domain into a number of subdomains in such a way that the summation of weights associated with the finite element nodes within each subdomain is approximately equal, and the communication between subdomains is minimised. With suitably chosen weights, it is expected that there will be approximately equal accuracy associated with each subdomain. This accuracy is maximised by allowing the partitioning to occur through areas of the domain that have relatively little flow activity, which, in this case, is characterised by the pointwise maximum Reynolds stresses. A Gaussian Process Regression (GPR) machine learning method is used to construct a set of local approximation functions (hypersurfaces) for each subdomain. Each local hypersurface represents not only the fluid dynamics over the subdomain it belongs to, but also the interactions of the flow dynamics with the surrounding subdomains. Thus, in this way, the surrounding subdomains may be viewed as providing boundary conditions for the current subdomain. We consider a specific example of turbulent air flow within an urban neighbourhood at a test site in London and demonstrate the effectiveness of the proposed DDNIROM. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The Reduced Order Model (ROM) approach is a very powerful numerical technique that can be extremely useful for real-time decision making and operational modelling as it enables predictions to be obtained in a time that is many orders of magnitude faster than conventional computational methods. Moreover, ROM solutions can be obtained within fractions of a second as opposed to days or weeks with conventional numerical methods. It is particularly popular in engineering and has been successfully applied in various areas of research such as fluid dynamics [1–6], data assim-

∗ Corresponding author at: ZCCE College of Engineering, Swansea University, Bay Campus, Fabian Way, Swansea, SA1 8EN, UK. E-mail addresses: [email protected], [email protected] (D. Xiao).

https://doi.org/10.1016/j.compfluid.2019.02.012 0045-7930/© 2019 Elsevier Ltd. All rights reserved.

ilation [7], fracture modelling [8], host-parasitoid system [9], air pollution [10], aerospace design [11], and haemodynamics [12]. ROMs can be classified in two groups based on the modification (or not) of the original source code; Non-Intrusive ROM (NIROM) refers to the ROM when the original source code is not modified, whilst intrusive ROM (IROM) refers to the ROM when the original source code is modified. NIROM has become particularly popular over the last decade as it has been shown to avoid instabilities and to overcome non-linear inefficiency issues encountered in IROM [13–18]. It also has the advantages over IROM of being easier to implement, modify and extend to use with complicated problems, as it does not rely on having access to the source code of the original model. Xiao et al. proposed a number of NIROMs [19–22] which have been applied within the field of fluid dynamics, fluid-structure interactions [23] and multiphase porous media flow [24].

16

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

Fig. 1. On-line NIROM calculation.

Fig. 2. Graph of the relationship between the maximum Reynolds stress error, τimax , for global NIROM and the maximum singular value λmax that is neglected after truncating i the SVD at a certain level. The graph also shows a fit to this data, given in Eq. (15).

However, despite its computational speed, reduced order modelling has certain difficulties. One of these is that it is hard to capture the details of a small region within a much larger computational domain e.g. capturing details of flow around a building in a town or a city. In order to overcome such deficiencies, methods which divide the domain into a number of subdomains have been applied to reduced order modelling. This idea of Domain Decom-

position (DD) goes back to the work presented by Schwarz [25] on overlapping domain decomposition. The idea was extended to nonoverlapping DD in the 1960s [26–28]. More recently there has been a lot of research carried out into combining DD methods with solvers, efficient preconditioners and parallel computing technologies, see for example [29–31]. For our purposes we exploit the basic idea of domain decomposition as our aim is to come up

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

17

Fig. 3. The plots show (a) a horizontal slice at a height of 15 m above the ground; (b) the surface mesh of the test site; (c) the mesh on a vertical slice through the centre-line of the tallest building and parallel to the streamwise direction.

with an effective way of splitting the domain. Pain et al. [32] also presented a domain decomposition method based on neural networks. Lucia et al. [33,34] firstly introduced the subdomain idea into ROM for tracking a moving shock wave, whilst other applications of the subdomain approach into ROM include the work of Baiges et al. [35], Kerfriden et al. [36], Amsallem et al. [37] and Chaturantabut [38]. Antil et al. [39,40] also used a subdomain idea to construct balanced truncation based ROMs. Xiao et al. [41] develop a Domain Decomposition NIROM (DDNIROM) which splits the computational domain into subdomains geometrically and calculates a local approximation to the governing equations within each subdomain. The DDNIROM is solved iteratively in order to pass information from one subdomain to another. DDNIROM can then be applied to large-scale computational problems for capturing complex flows in detail e.g. turbulent flows within the atmospheric boundary layer, flows with shocks, or flows at a city scale. DDNIROM has the advantage of allowing the construction of a set of local basis functions based on details of local forward flow solutions. This implies that the subdomains with less complicated physics can be modelled using fewer basis functions, whilst the others with more complicated flow patterns may be represented with a larger number of basis functions to enable capturing of the detailed flow and the turbulent eddies. However, none of these methods optimise the partitioning for the ROM according to the physics of the problem. Therefore, we extend our previous DDNIROM approach [41] to do just this. What distinguishes the current work from previous work on DD with ROM is that it carefully chooses the subdomains so that the accuracy of the ROM is improved. This is achieved through a weighting constraint which aims to achieve an equal accuracy in each subdomain as well as to minimise the dynamic activity between subdomains. The domain is decomposed so that the work load is balanced equally across each subdomain and the communication load (dynamic activity) between subdomains is minimised. Choosing weights that represent both the communication load and the load on each subdomain allows the subdomains to be more isolated (due to minimising the communication). The subdomains will not be independent (completely isolated) but perhaps more weakly linked than a random domain decomposition would produce. This careful choice of subdomains, and their boundaries, is the key to the success of the proposed method. This approach can be used for both intrusive and non-intrusive ROMs.

The main advantages of the DDNIROM approach can be summarised as follows: (a) it may be used for large-scale (km scale) problems with many degrees of freedom; (b) it is characterised by enhanced accuracy for the local hypersurface formations; (c) smaller Singular Value Decompositions (SVDs) within each subdomain as opposed to the larger SVDs required for the entire domain (SVDs are used to generate the basis functions); (d) Gaussian Process Regression (GPR), used in DDNIROM, being lower dimensional surface fits than those used in NIROM, may be more accurate because they surface fit in lower dimensions; (e) the potential ability to form the NIROMs for each subdomain independently by running models across a subdomain or group of subdomains independently. This makes possible the large-scale modelling of, for example, an entire city, even when the forward model is based on expensive computational methods. The structure of the present paper is as follows: Section 2 presents the governing equations of fluid flow problems. Section 3 describes the theory of the NIROMs. Section 4 derives the domain decomposition based NIROM (DDNIROM) approach. Section 5 demonstrates the performance of the DDNIROM for a specific turbulent flow test case: urban flows. Finally, in Section 6, the conclusions are presented. 2. Governing equations The three-dimensional continuity and non-hydrostatic Navier– Stokes equations which describe the conservation of mass and momentum of a fluid are given by

∇ · u = 0,

(1)

∂u + u · ∇ u = −∇ p + ∇ · τ , ∂t

(2)

where u ≡ (u, v, w)T is the velocity vector, t is the time, p = p˜/ρ0 is the normalised pressure ( p˜ being the pressure and ρ 0 being the constant reference density), and τ denotes the stress tensor. An anisotropic Smagorinsky model [42,43] is used to calculate the large eddy scale (LES) sub-grid scale viscosity which is included in the stress tensor. We apply filtering to the Navier–Stokes equations to model behaviour on the fine scale, for instance, fluctuations that occur on scales smaller than the grid scale.

18

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

Fig. 4. Plots of Reynolds stresses and the resulting domain decomposition. (a) and (c) show the Reynolds stresses on vertical and horizontal slices through the domain. In (a) the plane passes through the tallest building and is aligned with the streamwise direction. In (c) the plane lies 1 m above ground level. Plots (b) and (d) give an indication of some of the 32 different subdomains on a horizontal slice through the domain. Plot (b) zooms in on the central part.

3. Non-intrusive reduced order modelling This section describes the construction and use of NIROM for three dimensional urban flows using Proper Orthogonal Decomposition (POD) and the GPR machine learning method. Typically there are two stages linked to constructing a reduced order model: off-line (training) and on-line. The off-line stage, which is expected to be computationally expensive, involves finding the POD basis functions and approximating the governing equations in reduced space. The on-line stage is much less computationally intensive due to the reduced dimension of the model. During this stage the governing equations in reduced space are solved, which, in this case, amounts to mere function evaluations. The computational cost of these stages varies, depending on the application: we will give specific figures in the results section. NIROM can be deployed in two different ways, either to predict the flow at times beyond those used in the training of the model, or to predict the flow for parameters not used in the training. To explore the first case, we previously developed a NIROM for turbulent urban flows [44]. We ran the high-fidelity model until the flow statistics had reached a steady state and we showed that, when predicting beyond the time of the training period, the NIROM retained the same flow statistics as the high-fidelity model. In this paper we also run the DDNIROM beyond the training time period

and monitor the flow statistics. In the future we aim to develop parametrised NIROMs. In this section we first describe how POD is used to obtain basis functions from the high-fidelity LES model. We give a brief description of GPR and then explain how the NIROM is constructed by training a neural network using GPR with data from the highfidelity model. 3.1. Proper orthogonal decomposition POD, also known as principal component analysis (PCA), or empirical orthogonal functions (EOF) in oceanography, is one of the most widely used model reduction techniques. The aim of POD is to find a subspace that is able to approximate a given data set in an optimal least-squares sense, i.e. the method is based on finding a subspace from a full space such that the error of projecting from the full space onto the subspace is minimised [45]. In order to sample the high dimensional space we use the Method of Snapshots [46]. Snapshots are solutions of the discretised governing equations at certain time levels and are vectors of length N, where N is the number of nodes. If we have taken N s snapshots of the velocities, we construct snapshot matrices which have the following form



S = u1

u2 · · ·

uN s



(3)

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

19

Fig. 5. Plots showing the singular values versus number of POD basis functions with a logarithmic scale for: (a) global NIROM; (b) DDNIROM where the decomposition has been performed with a uniform nodal weighting, see Eq. (18); (c) DDNIROM a decomposition formed with weights based on the Reynolds stresses as defined in Eq. (16) (d) DDNIROM with the decomposition formed by weights based on an estimation of the largest singular value associated with the neglected modes in the SVD, see Eq. (15). The original NIROM result is shown in black on every plot and the other 32 subdomain results are shown in different colours.

where the discretised solutions ui are vectors containing values of the streamwise velocity at the nodes for the ith snapshot. An SVD is then applied to S which results in three matrices,

S = U V T

(4)

where U is an N by N matrix containing the left singular vectors, V is an N s by N s matrix containing the right singular vectors and  is an N by N s matrix which contains the singular values on the leading diagonal. The columns of matrix U contain the POD basis functions or modes. The singular values give an indication as to how important each basis function is, and the square of each singular value represents the amount of energy captured by the corresponding basis function, so if the singular value is low its corresponding basis function can be discarded. Based on this, a common criterion for choosing how many POD basis functions to include is as follows: for a tolerance η  1, find the smallest integer m such that

m

σi2  η, σ2 i=1 i i=1

N s

(5)

where σ i represents the ith singular value and m  N s (although the method works well when m  N s ). Considering the criterion, for example, if η takes the value 0.99, this means that 99% of the energy of the snapshots would be captured by the m leading POD basis functions. We find POD basis functions and singular values for each of the three velocity components in turn, because, in this work, the reduced order model will be formulated in terms of velocity alone.

As the problem is incompressible we can eliminate pressure from the NIROM formulation and pressure POD basis functions are not required. In this paper we obtain the snapshots of velocity by solving the discretised governing equations. After performing the SVD we have N s possible POD basis functions. After truncation, based on the energy criterion given in Eq. (5), we will have a set of m POD basis functions which we represent as {φ j }m . An approximation to the j=1 velocity can be expressed as a linear combination of the POD basis functions, m 

α jφ j

(6)

j=1

where α j represents the POD coefficients. 3.2. Gaussian process regression Gaussian Process Regression (GPR) (see [47]) uses a series of Gaussian-shaped basis functions to provide the surface representations used here to form the NIROM. We note that GPR produces smooth hypersurfaces which we believe regularises the solution to some extent, thereby increasing the accuracy of the results, see [48]. GPR often does not require as much training data as feed forward Neural Networks, and thus, because we have only 500 data points, we use GPR rather than feed forward Neural Networks. We adopt the GPR from within the open-source scikit-learn library [49] to form our DDNIROMs. Especially for small data sets, GPR can be susceptible to over-fitting, where the model learns the

20

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

Fig. 6. A comparison of velocity solutions between the high-fidelity model, NIROM with 48 POD bases at time instances t = 320 s (left panel) and t = 1400 s (right panel) on a horizontal plane at a height of 15 m above the ground.

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

21

Fig. 7. A comparison of velocity solutions between the high-fidelity model, NIROM with 48 POD bases and DDNIROM with 24 basis functions at time instances t = 320 s (left panel) and t = 1400 s (right panel) on a plane through the centre-line of the tallest building and parallel to the streamwise direction.

•

•

• •

run the high-fidelity LES model over the time period [0, T] and generate the snapshots matrix S; calculate a set of POD basis functions by performing the SVD of the snapshots matrix S; calculate the POD coefficients corresponding to each snapshot; train the GPR model using the POD coefficients of the snapshots, see Eqs. (8) and (9).

3.4. Using NIROM

Fig. 8. DDNIROM at time levels 2328 s, 2336 s, 2344 s, 2352 s on a plane through the centre-line of the tallest building and parallel to the streamwise direction.

During the online stage, the governing equations are approximated by the hypersurfaces { f j }m which are treated as response j=1 functions, allowing the NIROM solutions at the current time level to be predicted given those at a previous time level

α j (t + t ) = f j (α(t )) ∀ j ∈ {1, 2, . . . , m} .

(10)

training data well but predicts poorly for ‘unseen’ data. To reduce this effect, leave-p-out cross-validation is used, in which 80% of the data is used for training and the remaining 20% is used to judge the performance of the model (‘validation’). We allow the width of the Gaussian radial basis functions to be in the range [10−2 , 102 ] and the variance to be in the range [10−3 , 103 ]. Both quantities are optimised by the GPR algorithm.

We remark that when using the NIROM to predict, the time step, t, will necessarily coincide with the time between two successive snapshots. In this work we predict for times beyond those used in training the model, i.e. for times t > T. The procedure of using the NIROM is summarised in Fig. 1.

3.3. Training the non-intrusive reduced order model

4.1. Domain-decomposition-strategy

After obtaining the basis functions as described in Section 3.1, the neural network based on GPR is trained with the snapshot solutions from the LES simulations that have been projected onto the reduced space by the POD basis functions. The training process results in a set of hypersurfaces which will represent the governing equations, geometry, boundary conditions and model parameters in reduced space. The hypersurfaces are represented by a function fj for each POD coefficient, which maps the POD coefficients from one time level to the next, for example

This section describes a novel method that combines Domain Decomposition and NIROM (DDNIROM). It builds on a previous paper [41] in which DDNIROM was first developed. There, the decomposition was calculated simply by splitting the domain geometrically. The domain decomposition method used in the current paper calculates a balanced k-way partitioning that minimises the total communication volume [50] and distributes the computational work load as equally as possible over the subdomains. An important aspect is the overall accuracy of the DDNIROM which will be enhanced if the activity level between the domains is small measured by the maximum value in the Reynolds stress tensor at each point in the domain. In addition, although achieving an equal loading between subdomains may seem a secondary issue, as one always has the option of using variable numbers of POD basis functions within each subdomain, balancing the load within each subdomain may give us a suitable aim for both edge cut minimisation and load balancing. One can, of course, try to form iteratively an approximate weight set while finding the subdomains, by iteratively adjusting the weights until there is a balance in accuracy within each subdomain. We have not pursued this here as it could be complex and computationally expensive. Initially, weights are assigned to each vertex or node in the FE mesh and then the domain decomposition algorithm attempts to partition the domain into a number of subdomains, so that the

α kj = f j (αk−1 ) = f j (α1k−1 , α2k−1 , . . . , αmk−1 ),

∀k ∈ {1 , 2 , . . . , N s } . (7)

For each POD coefficient, the neural network is therefore trained with the following inputs and outputs

input:

output:

  αk−1 = α1k−1 , α2k−1 , . . . , αmk−1 α kj ,

(8)

(9)

for all k ∈ {1, 2, . . . , N s }. Through this training process we obtain a set of hypersurfaces (functions { f j }m ) which concludes the offj=1 line stage, the main steps of which can be summarised as follows:

4. Domain decomposition NIROM

22

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

Fig. 9. Subdomain number against number of nodes in each subdomain.

sum of the weights in the interface regions between subdomains is minimised and the sums of the weights in each of the subdomains are approximately equal. This has the effect of minimising the communication and activity between subdomains and balancing the computational load across the subdomains. We use the ParMETIS library [51] to decompose the computational domain. We quantify the communication load as follows. Consider a graph G = (V, E ) with a node set V and an edge set E. Let P be a vector of size |V| and let Pi store the subdomain number to which vertex i belongs. For each vertex i ∈ V let Ai be the number of subdomains other than Pi that the vertices adjacent to i belong. The total communication volume C of this partition is defined as:

C=



wi Ai ,

(11)

i∈V

see [51]. Eq. (11) corresponds to the total communication volume incurred by the partition because information associated with each interface vertex i needs to be sent to all of its Ai subdomains. In this work, we investigate three ways of calculating the vertex weights: uniform weights; weights based on the maximum value of the (nodal) Reynolds stresses, and weights based on an approximation of the singular value of the largest mode discarded in the truncation of the SVD. In the following, the x-component of the Reynolds Stresses defined as

u u =

1 T2 − T1

T2

(u − u )2 dt ,

(12)

T1

for a time interval of [T1 , T2 ], where u denotes the pointwise x-component (streamwise direction) of velocity and u denotes a time-averaged quantity. We first consider basing the nodal weights, wi , on the maximum Reynolds stress to attempt to minimise the activity at subdomain interfaces:

wi = max τi k l =: τimax . k,l

(13)

The computational domain is subsequently decomposed so that the amount information sent between the subdomains is minimised. We conjecture that as one minimises the activity between subdomains then one maximises the accuracy of the overall domain decomposition approach. We expect, further, that the domain decomposition will isolate subdomains in terms of the dynamics of the system. The second criterion (solutions in the subdomains have a similar accuracy), is also met by requiring that the sum of the weights be similar in each subdomain. Thus if the activity level within a subdomain is relatively low (as characterised by the

Fig. 10. Location of two points: red point 5 m above the ground level and the blue point is 10 m above the ground level. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Reynolds stresses) then this subdomain will contain more nodes than other subdomains with a greater activity level. Let us now consider basing the weights on the maximum singular value of the discarded modes of the SVD. We recall that the aim is to find a scalar field composed of the weights at each node, which can be used as a surrogate to approximate the communication between subdomains and the load within each subdomain. An equal load, activity or accuracy within each subdomain would be achieved if, for a given number of POD basis functions, the corresponding singular values that are neglected are approximately equal, i.e. the approximation to the ‘seen’ data within each subdomain is approximately the same. Using the original global NIROM model we can try to form a correlation between the maximum error in the Reynolds stresses τ err over the entire domain and the largest singular value λmax associated with the neglected POD modes. Looking at Fig. 2, we postulate that the relationship takes the following form:

  max τ err = α e(βλ ) − 1 .

(14)

Further we assume that (1) the maximum truncated singular value associated with each node can be calculated from the maximum value of the Reynolds stress tensor at each node, that is τ err = τ max ; and (2) that the sum of these singular values associated with each node should be balanced within each subdomain in order for there to be a balance of load or truncation of the singular values within each subdomain. Thus, setting the weights, at each node, equal to these singular values, at each node, is a suitable choice.

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

23

4.2. Domain decomposition NIROM This section describes how domain decomposition can be applied to NIROM. In the DDNIROM approach, the computational domain  is divided into S subdomains, d , d ∈ {1, 2, , S} and each d subdomain has local unknowns d ∈ RM , where Md is the number of nodes in subdomain d and M = M1 + · · · + Md + · · · + MS is the total number of nodes. In this method, the solutions at subdomain d are used to form a set of local POD basis functions d , the jth column of which is φdj . The local POD basis functions are calculated by performing an SVD on the nodal values of the snapshots from only this local subdomain. For each subdomain d , we construct a set of hypersurd faces { fid }m to represent the underlying dynamical system assoi=1 ciated with this subdomain and the surrounding subdomains over the reduced space (where md is the number of POD basis functions in subdomain d). Each hypersurface has the form of

αid,n = fid (αd,n−1 , αsd,n ),

(19)

where the vector αd,n−1 denotes the complete set of POD coefficients (i.e. the POD coefficients for all the solution variables required which is u, v and w in this case) at previous time level n − 1 for the subdomain d , αsd,n denotes the complete set of POD coefficients at time level n for the surrounding subdomains. The NIROM method is then used to construct a set of hypersurfaces for each POD coefficient of each subdomain. The procedure is described in Algorithm 1 (offline) and Algorithm 2 (online).

Fig. 11. The plots displayed above show the time series of the x-component of the velocity (m/s) at the two locations shown in Fig. 10 (the top plot corresponds to the higher (blue) point and the bottom plot to the lower (red) point) from the highfidelity model and from DDNIROM (in the predicting time interval) with 24 basis functions. The high fidelity model result is shown between 0 and 20 0 0 s and the DDNIROM is shown between 20 0 0 s and 40 0 0 s but on the same axis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Rearranging Eq. (14), and assuming that it holds at each node i, and using τ err = τ max gives

λmax = i

1

β



ln

τ

max i

α



+1 ,

(15)

in which α and β are two curve fitting scalars. We have found that α = 0.05 and β = 0.0019183 are suitable choices, see Fig. 2. We thus have three possible weights at each finite element node wi , that is:

wi = wi = or

τimax from Eq. (13 ), λmax from Eq. (15 ), i

wi = 1,

∀i.

(16) (17) (18)

It is worth noting that Eq. (17) is closer than Eq. (16) to the uniform weighting, i.e. it has less variation for a given Reynolds stress field. If we have other fields (other than just the velocity) then we might form similar expressions to the above based on the standard deviation of these fields rather than Reynolds stresses, say. We are then presented with choices on how we balance the subdomains based on two or more weights at each node. A maximum criteria would be one option, for example.

Algorithm 1: Offline: constructing a set of hypersurfaces for DDNIROM. (1) Generate Ns snapshots over the time period [0, T ] by solving the governing equations (2) Divide the computational domain  into S subdomains (3) Generate POD basis functions and coefficients for d = 1 to S do (a) Generate the POD basis functions d by performing an SVD on the snapshots matrix of subdomain d; (b) Calculate a set of POD coefficients (α1d,n , α2d,n , · · · , α d,n ) md for each subdomain by projecting the snapshots into a reduced space (md is the number of POD coefficients in subdomain d); endfor (4) Obtain a set of hypersurfaces for d = 1 to S do for i = 1 to md do 

(i) Set the values yid = αid,1 , αid,2 , · · · , αi

d,Nt

for the

ith POD coefficient (ii) Obtain the hypersurface fid , for subdomain d and POD coefficient i, by training the GPR network; endfor endfor

Algorithm 1 presents the offline computational procedure of how to construct a set of local hypersurfaces for each subdomain while Algorithm 2 describes the online computation of DDNIROM where interactions between a subdomain and its surrounding subdomains are taken into account. At every time level, information about the incompressibility must be communicated to all the subdomains,

24

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27 Table 1 L2 Error norms for NIROM (48 basis functions) and DDNIROM (24 basis functions) calculated for the streamwise velocity at 320 s using the high-fidelity model as a reference solution. (See Eq. (28) for the definition of the error norms.).

Algorithm 2: Online DDNIROM calculation for the fluid problem. for n = 1 to Nt do  time loop  calculate the POD coefficients at the current time step: for it = 1 to Niteration do  iteration loop

t = 320 s

t = 1400 s

eNIROM L2

24.78%

23.61%

eDDNIROM L2

10.98%

15.23%

for d = 1 to S do  subdomain loop if it == 1 Initialise POD coefficients for subdomain d (αd,0 ) for i = 1 to md do  POD coefficient loop (i) Evaluate the hypersurface fid by using the complete set of POD coefficients αd,n−1 and, if available, αsd,n , which gives the POD coefficient αid,n at the current time level n:

αid,n = fid (αd,n−1 , αsd,n )

(20 )

(ii) The steps inside this loop should be carried out for the POD coefficients of every variable, i.e. all the velocity components endfor endfor When it > 1 check convergence: (it )

errord,n =

||(it ) αd,n − (it−1) αd,n ||∞ < tol ||(it ) αd,n ||∞

(21 )

 project the POD coefficients to the full space for d = 1 to S do  Calculate the solution ud,n for each subdomain d for each time step n by projecting α d,n onto the full space. j

ud,n = u¯ d +

md 

α d,n φdj j

(22 )

In this work we use just one iteration and the results from the DDNIROM were found to be close to the full model (see Fig. 6 and Table 1). The reason why we found a single iteration was adequate, for this example, was that these subdomains have been defined so that there is relatively little activity between subdomains. This enables the iterative domain decomposition method to converge quickly. In addition, even a single iteration has a level of associated implicitness, that is, the solution is updated by using information from neighbouring subdomains as one sweeps across the subdomains. 5. Numerical example: modelling air flow around london south bank university This section demonstrates the ability of DDNIROM using a test case in the London South Bank University (LSBU) area, shown in Fig. 3. The computational domain has a size of [0, 2041] × [0, 2288] × [0, 250] (metres). A synthetic incoming-eddy method is used at the inlet [52]. The mean velocity profile is prescribed as in Eq. (25)

z   (u, v, w ) = 0.97561 ln , 0, 0 0.01

where z denotes the height (in m). The Reynolds stresses Re (Eq. (26)) and the length-scale L (Eq. (27)) are prescribed constant and equal to 0.8 and 100 m respectively, for the diagonal components, and zero elsewhere.



Re =

j=1



endfor endfor

L=

therefore we chose an implicit coupling resulting in a recurrent neural network. In Algorithm 2, the inputs of the function fid include POD coefficients from the previous time level αd,n−1 (known) and the current time level from surrounding subdomains αsd,n . The latter might be unknown, in which case, we take the value from the previous time level. In order to solve the implicit system we introduce an iteration method to obtain αd,n at current time level n. The iteration loop (for it = 1 to Niteration ) ensures that the hypersurface incorporates the flow dynamics over the subdomain d, but also, includes the flow interaction between the subdomain and its neighbouring subdomains. Convergence could be checked by using the following criterion (it )

errord,n =

||(it ) αd,n − (it−1) αd,n ||∞ < tol ||(it ) αd,n ||∞

(23)

where tol is a tolerance provided by the user and the infinity norm is defined as

||x||∞ = max x j . j

(24)

(25)

0.8 0 0

100 0 0



0 0.8 0

0 0 , 0.8

0 100 0

0 0 . 100

(26)



(27)

Zero velocity is prescribed on the wall boundaries and bottom. A perfect slip condition is specified on the vertical lateral boundaries and zero stress conditions is set to be p = 0 at the outlet boundary. In this case, the streamwise direction corresponds to a westerly wind direction. The high-fidelity solutions were obtained by running in parallel the open-source finite element fluid dynamics code called Fluidity [53,54]. Fluidity has been parallelised using standard domain decomposition techniques, see [55], where the computational domain is partitioned and each processor is assigned a partition. The Message Passing Interface library is used to communicate between processors. The velocity equations were solved with the Generalised Minimal Residual method and a symmetric successive over-relaxation preconditioner, and the pressure equations were solved with the conjugate gradient method and the BoomerAMG preconditioner. More information about these solvers and preconditioners can be found in the PETSc Users’ Manual [55] and technical report [56]. To initialise the problem, the LES simulation was performed for one hour (in real time). The mesh was adapted every three time steps, specifying a maximum edge length of 50 m and a minimum

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

25

Fig. 12. Probability density functions of the x-component of velocity of the air at 5 m and 10 m above the ground level and at the two points shown in Fig. 10. Left are the high fidelity model results and right the DDNIROM results with 24 POD basis functions in each subdomain. The horizontal axis is the interval of velocity (m/s) in the x-direction or streamwise direction.

edge length of 0.3 m. This time interval (one hour real time) is sufficient for the flow statistics to reach a quasi-steady state. From this point onwards, the mesh is fixed, as shown in Fig. 3. After the initialisation stage, the simulation is continued for 20 0 0 s with a time-step size of 4/3 s on the fixed unstructured mesh with 767,559 nodes. During the 20 0 0 s, 50 0 snapshots were taken every 4 s from the high-fidelity LES simulation results. The DDNIROM was trained with these 500 snapshots. Fig. 4 shows the x component of the Reynolds stress (u u ) on a vertical plane through the largest building in the domain and this plane is aligned with the streamwise direction. In addition, in the same figure, we show the same Reynolds stresses at a height of 1m above the ground and viewed from underneath the domain. We also show the some of the subdomains viewed from underneath the domain at this height. What should be noticed is that the subdomain boundaries tend to avoid the areas with large Reynolds stresses. These results and, unless stated, all that follow, use the maximum Reynolds stresses to determine the nodal weights in the decomposition. The singular values associated with the SVD of the global NIROM are shown in Fig. 5 and are slow to decrease after the

30th singular value. Also shown in this figure are the singular values associated with the subdomains following the DDNIROM approaches, with vertex weights determined by Eqs. (16)–(18). Notice that the singular values associated with the subdomains are always less than the global NIROM, but for the uniform weight case wi = 1 some of the subdomains have singular values close to the global NIROM. This suggests that in parts of the domain the solution will not be accurate for this DDNIROM. However, if one uses the Reynolds stresses or approximate singular values as weighting parameters (Fig. 5(c) and (d)) then the singular values associated with all the subdomains are substantially lower than those associated with the global NIROM with the weights defined by the singular values, Eq. (17), being slightly better (that is, decaying more rapidly). This is the reason why the resulting DDNIROMs are more accurate, see Figs. 6 and 7 obtained with only 24 basis functions in each subdomain. Notice in these images that the DDNIROM with 24 basis functions is considerably more accurate than the global NIROM with 48 basis functions at the time levels shown and, in fact, at all times. The L2 error norm of the streamwise velocity is presented in Table 1 for both NIROM and DDNIROM, in which the high-fidelity model is taken as the reference solution. The norm is

26

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27 Table 2 Comparison of the CPU cost (in seconds) required to solve the high fidelity model and DDNIROM for 4 s in real time. For the DDNIROM, ‘solution’ corresponds to step (a) in the flow chart, Fig. 1, and ‘projection’ corresponds to step (b).

Assembly Solution Projection Total

High-fidelity model

DDNIROM

}1555

n/a 0.128 0.003 0.131

n/a 1555

RAM. As the DDNIROM has so few degrees of freedom there would be no gain from running this model in parallel. The total amount of time required to generate the 500 snapshots for this example was 9 days approximately. 6. Conclusions

defined as follows:

eNIROM (t ) = L2

||uHFM (t ) − uNIROM (t )||2 ||uHFM (t )||2

(28)

where uNIROM (t) represents the streamwise nodal velocities from NIROM at time t and uHFM (t) represents the streamwise nodal velocities from the high-fidelity model at time t. The error norm for DDNIROM is similarly defined. In Fig. 7 we compare the flow speed of the high-fidelity model, the NIROM with 48 POD bases and the DDNIROM with 24 basis functions at time instances t = 320 s (left panel) and t = 1400 s (right panel). The results are shown on a plane through the centreline of the tallest building and parallel to the streamwise direction. This figure shows the superior performance of DDNIROM over NIROM with 48 basis functions as well as the ability of DDNIROM to reproduce closely the high-fidelity model result. In Fig. 8 we see a similar plot but of DDNIROM at a number of consecutive time levels. In this plot can we can see that the DDNIROM is able to capture eddies moving downstream of the tall building. In Fig. 9 we show the number of nodes within each subdomain for the three methods of determining weights for the decomposition, given by Eqs. (16)–(18). In this figure the subdomains are ordered in increasing number of nodes belonging to them. Notice that the method with equal weights Eq. (18) has approximately equal number of nodes within each subdomain. This is approximate because the ParMETIS partitioning method also tries to minimise the communication between subdomains as defined by Section 4 and thus balancing the weight sum in each subdomain is only approximate. Also notice that the weights based on the maximum stress, Eq. (16), result in the biggest variation of the node numbers in each subdomain simply because it has the biggest variation between Eqs. (16) and (17). This increased variation can be seen in Fig. 2. In Fig. 11 we show the time series of two points in the domain (near the busy intersection, see Fig. 10). Here we show on the same graph the first 20 0 0 s obtained from the LES model and results from DDNIROM from the range [20 0 0,40 0 0]. Visually the LES and DDNIROM models seem to reproduce similar frequencies and amplitudes of the x-component of velocity (streamwise component). Moreover, they also produce similar probability density functions, see Fig. 12. This helps to show that the DDNIROM model is an accurate model in itself with similar dynamics to the high fidelity model. Table 2 shows the CPU cost required by both the high-fidelity model and DDNIROM to solve for 4 s in real time (which corresponds to three time steps for the high-fidelity model and one time step for the DDNIROM). The high-fidelity model solves for just over 3 million degrees of freedom, whereas the DDNIROM has just 72 degrees of freedom in each subdomain. It is worth noting that the CPU time required to solve the DDNIROM over this time interval is only 0.128 s, whereas the high-fidelity model requires 1555 s running in parallel on a workstation with 10 Intel®Xeon®X5680 CPU processors (each core has a frequency of 3.3 GHz) and 512GB

In this paper we develop a method of decomposing the computational domain into subdomains in order to enhance the performance of the Domain Decomposition Non-Intrusive Reduced Order Model (DDNIROM). Within each subdomain we use the Gaussian Process Region (GPR) (a machine learning method) to obtain a hypersurface that approximates the discretised governing equations. Each local hypersurface represents, not only the fluid dynamics over the subdomain it belongs to, but also the interactions with the surrounding subdomains. This implicit coupling between the subdomains provides the global coupling necessary to enforce incompressibility and is a means of providing boundary conditions for each subdomain. As suggested, choosing the subdomains is the key to the good performance of the method and we have investigated three different choices of weighting (of the finite element nodes) for the domain decomposition algorithm. If we use a node weight that is based on the maximum Reynolds stresses or the largest singular value of the discarded POD modes, then the subdomain approach worked well, because, by design, there is relatively little activity between subdomains. However, if one uses a uniform weight then the accuracy of the approach was not as good, and only marginally better than global NIROM (DDNIROM with one subdomain). In addition, it may not be so important to balance the accuracy associated with a fixed number of POD basis functions in each subdomain. The reason for this is because we can always choose a different number of basis functions associated with each subdomain as demonstrated previously in [41] in order to balance the accuracy within each subdomain. The major advantage of using domain decomposition methods is that they: (a) may be used for large problems with lots of degrees of freedom; (b) are more accurate for the local hypersurface formation; (c) allow the use of smaller Singular Valued Decompositions within each subdomain when compared to SVDs across the entire domain; (d) lead to lower dimensional surface fits in the GPR than used in NIROM, which may result in a higher accuracy because they surface fit in lower dimensional spaces; (e) result in an ability to form the NIROMs independently by running models across each region independently. This may make possible large-scale modelling, of, say an entire city, even based on detailed LES turbulent flow models. These advantages have been demonstrated within this paper. In addition, it is shown that even using relatively small number of POD basis functions within each subdomain then it is possible to form a high fidelity DDNIROM: only 24 basis functions were used in each subdomain in this work. The reason for this is the choice of the subdomains, in particular, minimising the activity between subdomains. This was demonstrated in the complex problem of turbulent flows around buildings. This is a fairly typical turbulent flow problem and thus we expect similar performance in other turbulent flow problems which will be the subject of future investigations. Acknowledgments The authors would like to thank the reviewers for their helpful comments which have improved the manuscript. The authors are grateful for the support of the following: Managing Air for Green Inner Cities (MAGIC) (EP/N010221/1); the EPSRC multi-phase flow

D. Xiao, C.E. Heaney and F. Fang et al. / Computers and Fluids 182 (2019) 15–27

programme grant MEMPHIS (EP/K003976/1); The IMPACT (Innovative Materials, Processing and Numerical Technologies) at Swansea University; the Innovate UK Smart-GeoWells consortium; Imperial College High Performance Computing Service; funding from the European Union Seventh Frame work Programme (FP7/20072013) under grant agreement No.603663 for the research project PEARL (Preparing for Extreme And Rare events in coastal regions) and the EPSRC MUFFINS project (MUltiphase Flow-induced Fluid-flexible structure InteractioN in Subsea applications EP/P033148/1). References ´ [1] Noack BR, Stankiewicz W, Morzynski M, Schmid P. Recursive dynamic mode decomposition of transient and post-transient wake flows. J Fluid Mech 2016;809:843–72. [2] Xie X, Wells D, Wang Z, Iliescu T. Approximate deconvolution reduced order modeling. Comput Methods Appl Mech Eng 2017;313:512–34. [3] Lorenzi S, Cammi A, Luzzi L, Rozza G. POD-Galerkin method for finite volume approximation of Navier–Stokes and RANS equations. Comput Methods Appl Mech Eng 2016;311:151–79. [4] Wang Y, Navon IM, Wang X, Cheng Y. 2D Burgers equation with large Reynolds number using POD/DEIM and calibration. Int J Numer Methods Fluids 2016;82(12):909–31. [5] Noack BR, Morzynski M, Tadmor G. Reduced-order modelling for flow control, 528. Springer; 2011. [6] Bergmann M, Cordier L. Optimal control of the cylinder wake in the laminar regime by trust-region methods and POD reduced-order models. J Comput Phys 2008;227(16):7813–40. [7] Dimitriu G, Apreutesei N, S¸ tefa˘ nescu R. Numerical simulations with data assimilation using an adaptive POD procedure. In: Large-Scale Scientific Computing. Springer; 2010. p. 165–72. [8] Oliver J, Caicedo M, Huespe A, Hernández J, Roubin E. Reduced order modeling strategies for computational multiscale fracture. Comput Methods Appl Mech Eng 2017;313:560–95. [9] Dimitriu G, Stefanescu R, Navon IM. POD-DEIM approach on dimension reduction of a multi-species host-parasitoid system. Ann Acad Rom Sci Ser Math Appl 2015;7(1):173–88. [10] Fang F, Zhang T, Pavlidis D, Pain C, Buchan AG, Navon IM. Reduced order modelling of an unstructured mesh air pollution model and application in 2D/3D urban street canyons. Atmos Environ 2014;96:96–106. [11] Manzoni A, Salmoiraghi F, Heltai L. Reduced Basis isogeometric methods (RB-IGA) for the real-time simulation of potential flows about parametrized NACA airfoils. Comput Methods Appl Mech Eng 2015;284:1147–80. [12] Ballarin F, Faggiano E, Ippolito S, Manzoni A, Quarteroni A, Rozza G, et al. Fast simulations of patient-specific haemodynamics of coronary artery bypass grafts based on a POD–Galerkin method and a vascular shape parametrization. J Comput Phys 2016;315:609–28. [13] Schlegel M, Noack B. On long-term boundedness of Galerkin models. J Fluid Mech 2015;765:325–52. doi:10.1017/jfm.2014.736. [14] Osth J, Noack BR, Krajnovi S, Barros D, Bore J. On the need for a nonlinear subscale turbulence term in POD models as exemplified for a high-Reynoldsnumber flow over an Ahmed body. J Fluid Mech 2014;747:518–44. doi:10.1017/ jfm.2014.168. [15] Franca LP, Frey SL. Stabilized finite element methods: II. The incompressible Navier–Stokes equations. Comput Methods Appl Mech Eng 1992;99(2):209–33. [16] Chaturantabut S, Sorensen DC. Nonlinear model reduction via discrete empirical interpolation. SIAM J Sci Comput 2010;32:2737–64. [17] Xiao D, Fang F, Du J, Pain CC, Navon IM, Buchan AG, et al. Non-linear Petrov–Galerkin methods for reduced order modelling of the Navier–Stokes equations using a mixed finite element pair. Comput Methods Appl Mech Eng 2013;255:147–57. [18] Xiao D, Fang F, Buchan AG, Pain C, Navon IM, Du J, et al. Non-linear model reduction for the Navier–Stokes equations using residual DEIM method. J Comput Phys 2014;263:1–18. [19] Xiao D, Fang F, Buchan AG, Pain C, Navon IM, Muggeridge A. Non-intrusive reduced order modelling of the Navier–Stokes equations. Comput Methods Appl Mech Eng 2015;293:522–41. [20] Xiao D, Fang F, Pain C, Navon IM. A parameterized non-intrusive reduced order model and error analysis for general time-dependent nonlinear partial differential equations and its applications. Comput Methods Appl Mech Eng 2017;317:868–89. [21] Wang Z, Xiao D, Fang F, Govindan R, Pain C, Guo Y. Model identification of reduced order fluid dynamics systems using deep learning. Int J Numer Methods Fluids 2018;86(4):255–68. doi:10.1002/fld.4416. [22] Xiao D, Fang F, Pain C, Hu G. Non-intrusive reduced order modelling of the navier-stokes equations based on RBF interpolation. Int J Numer Methods Fluids 2015;79(11):580–95. doi:10.1002/fld.4066.

27

[23] Xiao D, Yang P, Fang F, Xiang J, Pain C, Navon IM. Non-intrusive reduced order modeling of fluid-structure interactions. Comput Methods Appl Mech Eng 2016;303:35–54. [24] Xiao D, Lin Z, Fang F, Pain C, Navon IM, Salinas P, et al. Non-intrusive reduced-order modeling for multiphase porous media flows using smolyak sparse grids. Int J Numer Methods Fluids 2017;83(2):205–19. [25] Schwarz H. Ueber einige abbildungsaufgaben. Journal für die reine und angewandte Mathematik 1869;70:105–20. [26] Przemieniecki J. Matrix structural analysis of substructures. AIAA 1963;1(1):138–47. [27] Schur I. Gesammelte Abhandlungen. Band II.; 1973. [28] Schur I. Vorlesungen über invariantentheorie; 1968. [29] Meurant G. Domain decomposition methods for partial differential equations on parallel computers. Int J Supercomput Appl 1988;2(4):5–12. [30] Quarteroni A, Valli A. Domain decomposition methods for partial differential equations. Oxford Science Publications; 1999. [31] Saad Y. Iterative methods for sparse linear systems. SIAM; 2003. [32] Pain C, Goddard A. A neural network graph partitioning procedure for grid based domain decomposition. Int J Numer Methods Eng 1999;44:593–613. [33] Lucia D, King P, Beran P. Reduced order modeling of a two-dimensional flow with moving shocks. Comput Fluids 2003;32(7):917–38. [34] Lucia D, King P, Beran P. Domain decomposition for reduced-order modeling of a flow with moving shocks. AIAA J 2002;40(11):2360–2. [35] Baiges J, Codina R, Idelsohn S. A domain decomposition strategy for reduced order models. Application to the incompressible Navier–Stokes equations. Comput Methods Appl Mech Eng 2013;267:23–42. [36] Kerfriden P, Goury O, Rabczuk T, Bordas S. A partitioned model order reduction approach to rationalise computational expenses in nonlinear fracture mechanics. Comput Methods Appl Mech Eng 2013;256:169–88. [37] Amsallem D, Zahr M, Farhat C. Nonlinear model order reduction based on local reduced-order bases. Int J Numer Methods Eng 2012;92(10):891–916. [38] Chaturantabut S. Temporal localized nonlinear model reduction with a priori error estimate. Appl Numer Math 2017;119:225–38. [39] Antil H, Heinkenschloss M, Hoppe R, Sorensen D. Domain decomposition and model reduction for the numerical solution of PDE constrained optimization problems with localized optimization variables. Comput Vis Sci 2010;13(6):249–64. [40] Antil H, Heinkenschloss M, Hoppe R. Domain decomposition and balanced truncation model reduction for shape optimization of the Stokes system. Optim Methods Softw 2011;26(4–5):643–69. [41] Xiao D, Fang F, Heaney C, Navon IM, Pain C. Domain decomposition for the non-intrusive reduced order modelling of fluid flow. Comput Methods Appl Mech Eng 2019. submitted. [42] Aristodemou E, Bentham T, Pain C, Colvile R, Robins A, ApSimon H. A comparison of mesh-adaptive LES with wind tunnel data for flow past buildings: Mean flows and velocity fluctuations. Atmos Environ 2009;43(39):6238–53. [43] Song J, Fan S, Lin W, Mottet L, Woodward H, Wykes MD, et al. Natural ventilation in cities: the implications of fluid mechanics. Build Res Inf 2018;46(8):809–28. doi:10.1080/09613218.2018.1468158. [44] Xiao D, Heaney C, Mottet L, Fang F, Lin W, Navon IM, et al. A reduced order model for turbulent urban flows using machine learning. Build Environ 2019;148:323–37. [45] Pinnau R. Model reduction via proper orthogonal decomposition. Springer Berlin Heidelberg; 2008. p. 95–109. [46] Sirovich L. Turbulence and the dynamics of coherent structures, part III: dynamics and scaling. Q Appl Math 1987;XLV:583–90. [47] Rasmussen C, Williams C. Gaussian processes for machine learning. The MIT Press; 2006. [48] Rasmussen C. Gaussian processes in machine learning. In: Advanced lectures on machine learning. Springer; 2004. p. 63–71. [49] Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: machine learning in python. J Mach Learn Res 2011;12:2825–30. [50] Karypis G, Kumar V. Multilevel k-way partitioning scheme for irregular graphs. J Parallel Distrib Comput 1998;48(1):96–129. [51] Karypis G, Schloegel K, Kumar V. ParMetis: Parallel Graph Partitioning and Sparse Matrix Ordering Library; 1997. [52] Pavlidis D, Gorman GJ, Gomes JLMA, Pain CC, ApSimon H. Synthetic-eddy method for urban atmospheric flow modelling. Boundary Layer Meteorol 2010;136:285–99. [53] Pain C, Piggott M, Goddard A, Fang F, Gorman G, Marshall D, et al. Three-dimensional unstructured mesh ocean modelling. Ocean Modell 2005;10(1–2):5–33. [54] Fluidity Version 4.1.12; 2015. 10.6084/m9.figshare.1387713. [55] Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, et al. PETSc Users Manual Tech. Rep. ANL-95/11 - Revision 3.10. Argonne National Laboratory; 2018. [56] Balay S, Gropp WD, McInnes LC, Smith BF. Efficient management of parallelism in object oriented numerical software libraries. In: Arge E, Bruaset AM, Langtangen HP, editors. Modern Software Tools in Scientific Computing. Birkhauser Press; 1997. p. 163–202.