Michele Quattromini
email: michele.quattromini@poliba.it

Michele Alessandro Bucci

Stefania Cherubini
email: stefania.cherubini@poliba.it

Onofrio Semeraro
email: onofrio.semeraro@universite-paris-saclay.fr

OPERATOR LEARNING OF RANS EQUATIONS: A GRAPH NEURAL NETWORK CLOSURE MODEL A PREPRINT

The spread of machine learning (ML) techniques in combination with the availability of high-quality experimental and numerical data boosted in recent years numerous applications in fluid mechanics. Among those, examples of closure models for turbulent flows or data-assimilation based on neural networks (NN) are already numerous. However, it is well known that these techniques are prone to over-fit and necessitate an exceedingly large amount of data, unless enforcing physical constraints. We address those limitations by applying graph neural networks (GNN). This architecture is characterized by a net of nodes that can be easily adapted to unstructured meshes. Moreover, it is known GNN can show remarkable generalization capabilities as compared to standard network models. Here, we demonstrate the use of GNN by interfacing them with a finite elements (FEM) solver for the supervised learning of Reynolds-averaged Navier-Stokes equations. We consider as a test-case the flow past bluff bodies; we train the model using the wake past a cylinder, at different Reynolds numbers 40 ≤ Re ≤ 150 and resolved on different grids. The GNN model successfully predicts the Reynolds stress also in unseen cases characterized by different geometric configurations and Reynolds numbers. Interestingly, a small data-set is used for achieving these performances, thus suggesting the applications of GNN in replacement of less flexible techniques such as convolutional NN.

Introduction

Machine learning (ML) applications are spreading across the scientific community in the most diverse fields, leading to the development of novel techniques and algorithms encompassed within the emerging field of scientific machine learning. This breakthrough is fostered by the massive augmentation of computer-aided activities, the constant increase of high performance computational power, the recent developments in the field of deep learning and the large availability of data [START_REF] Goodfellow | Deep learning[END_REF].

In a very schematic way, ML techniques enable to identify mappings between observables of a system (inputs) and quantities of interest (outputs) we aim to predict by leveraging data; when these analysed data are governed by deterministic or statistical laws, in principle, these mappings correspond to approximating models. Due to these potentialities, the number of contributions combining ML and fluid mechanics has been constantly growing in a large variety of applications, ranging from modelling of complex behaviour to reinforcement learning [START_REF] Brunton | Machine learning for fluid mechanics[END_REF]. These applications found a natural ground into modelling for turbulent flows. In [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF], a broad overview can be found on the different levels of approximation, together with a critical take on the limitations of the approach; in the review by [START_REF] Ghattas | Learning physics-based models from data: perspectives from inverse problems and model reduction[END_REF], the focus is on inverse modelling and model reduction. Among the numerous authors that addressed this problem, [START_REF] Ling | Evaluation of machine learning algorithms for prediction of regions of high reynolds averaged navier stokes uncertainty[END_REF] applied classification methods for identifying regions of uncertainties where the closure term of the Reynolds-averaged Navier-Stokes (RANS) might fail; [START_REF] Ströfer | End-to-end differentiable learning of turbulence models from indirect observations[END_REF] combined data assimilation with neural networks (NN) modelling of the Reynolds stress using limited observation. Other approaches leverage baseline models such as the the Spalart-Allmaras closure [START_REF] Singh | Using field inversion to quantify functional errors in turbulence closures[END_REF], physics-informed NN (PINNs) [START_REF] Eivazi | Physics-informed neural networks for solving reynolds-averaged navier-stokes equations[END_REF], regression methods [START_REF] Schmelzer | Machine learning of algebraic stress models using deterministic symbolic regression[END_REF] or decision trees [START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF].

In principle, we could identify universal closure models directly from data, having at disposal an infinite amount of them; in practice, we often deal with a limited amount of data or with few measurements available in the domain of interest. This impacts on the use of methods such as NN, where the expressivity and generic structure make them suitable for a large class of models, yet prone to generalisation problems. In this sense, ML models risk to be representative solely of the datasets included in the training process; thus, it becomes compelling given the available data to circumvent these problems by inputting well selected data during the training or providing prior knowledge through modelling [START_REF] Bucci | Leveraging the structure of dynamical systems for data-driven modeling[END_REF], [START_REF] Shukla | Scalable algorithms for physics-informed neural and graph networks[END_REF]. Alongside those limitations, when unstructured meshes are considered, standard techniques such as convolutional NN (CNN) might be limited or would require an interpolation step.

In this paper, we address some of these limitations by introducing graph neural networks (GNN) [START_REF] Hamilton | Graph representation learning[END_REF] in combination with finite elements method (FEM) for the simulations of unsteady flows. GNNs have a strong potential due to their peculiarities: these architectures are characterized by complex multi-connected nets of nodes that can be provided as input and easily adapted to unstructured meshes. Indeed, the convolution in a GNN is performed by aggregating information from neighbouring nodes, thus overcoming the limitations imposed by the geometry in contrast with CNN. Moreover, besides being differentiable, it is known GNN models show remarkable generalization capabilities as compared to standard network models [START_REF] Sanchez-Gonzalez | Graph networks as learnable physics engines for inference and control[END_REF]. Finally, the possibility of directly targeting the learning of the operator via the GNN's discrete stencils has been subjected to research (see [START_REF] Shukla | Scalable algorithms for physics-informed neural and graph networks[END_REF]). Among these works, we take inspiration from the contribution by [START_REF] Donon | Deep Statistical Solvers[END_REF], where the authors combined the GNN-based architecture -for incorporating permutation and translation invariance -with the statistical solver problem, and proof that this method has some universal approximation properties. Thus, GNN can be combined with unsupervised learning to perform operator learning for the solution of partial differential equations (PDE), while retaining a great flexibility with respect to the domain of computation. Following this philosophy, but using supervised learning, we consider the dynamics of the wake past a cylinder, at low Reynolds numbers 40 ≤ Re ≤ 150 and infer a closure model of the associated RANS equations, which generalizes on unseen Reynolds numbers and geometries. More in detail, a GNN model is trained on a dataset composed by the meanflows and the Reynolds stress of the wake past a cylinder, computed for O(10) cases at different Re on meshes characterized by different levels of refinement. We show that such a model reproduces with satisfactory accuracy the Reynolds stress of unseen cases. Moreover, leveraging the flexibility of the GNN, the model is capable at predicting -although with a more limited accuracy -the Reynolds stress on unseen geometries and flow conditions, thus demonstrating the potentiality of operator learning.

The remainder of the article is organised as follows: in Sec. 2 we introduce the flow cases and parameters; Sec. 3 describes the data-driven strategy; results are illustrated in Sec. 4 and discussed in Sec. 5.

Governing equations and numerical simulations

We consider incompressible two-dimensional (2D) fluid flows developing past bluff bodies. Here, we will focus on time-averaged quantities and the second order-statistics, the Reynolds stress. To this end, we introduce the Reynolds decomposition

u(x,t) = u(x) + u (x,t), (1)
where u = (u, v) is the meanflow and u the fluctuation field. Formally, any unsteady flow can be described through this decomposition, whether we consider a laminar or turbulent case [START_REF] Foures | A data-assimilation method for reynolds-averaged navierstokes-driven mean flow reconstruction[END_REF]. Plugging Eq. ( 1) in the Navier-Stokes equations and averaging we get the system of equations

u • ∇u + ∇p - 1 Re ∇ 2 u = f (2a) ∇ • u = 0, ( 2b
)
where p is the average pressure field. The Reynolds number is defined as Re = U ∞ D/ν, with the reference velocity chosen as the free stream velocity U ∞ at the inlet, D the reference length for the analysed flow and ν the kinematic viscosity. Mathematically, the resulting system provides the Reynolds-averaged Navier-Stokes Equations (RANS) and the forcing f is the closure term or Reynolds stress. This term can be modelled or -when data are availabledirectly computed as f = -u • ∇u .

(3) The forcing term f can be assumed as an unknown variable of a data assimilation scheme; if the RANS formulation is adopted as a baseline model, the identification of the forcing becomes a control parameter of the optimization process [START_REF] Foures | A data-assimilation method for reynolds-averaged navierstokes-driven mean flow reconstruction[END_REF].

Case Study design: bluff body flows

As a reference case of a bluff body flow, we consider the unsteady wake developing past a cylinder at Reynolds numbers 40 ≤ Re ≤ 150. This flow is commonly found as a benchmark and well documented: it exhibits stable behaviour up to a critical Reynolds number of Re c ∼ = 46.7, when a supercritical Hopf bifurcation occurs [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. At higher Reynolds numbers, the baseflow becomes an unstable solution, while the unsteady flow develops into a limit cycle behaviour known as von Karman street. We will focus on values of Re exhibiting this dynamic. In Fig. 1, the geometry of the test case is shown. The diameter of the cylinder is kept as the characteristic size of the analysed case, D = 1. Based on this dimension, the domain extension is L x = 30 in the streamwise and L y = 10 in the transverse direction; in the basic configuration, the cylinder is placed at (0, 0), at a distance ∆x = 11 from the inlet. The flow evolves from the inlet, where a uniform velocity u = (1, 0) T is imposed. At the outlet, the pressure is null, while on the top and bottom boundaries we set ∂ y u = 0 and v = 0. No slip conditions are considered around the obstacle. All the numerical simulations are initiated with null flow fields at t = 0 and performed using the FEM approach, supported by the FEniCS library [START_REF] Alnaes | The fenics project version 1.5[END_REF]; more details are provided in the supplementary appendix. In Fig. 2a, the meanflow along with the isolines of the vorticity ω are shown for Re = 150. Statistics are computed on-the-fly during the simulation. The final time of simulations is determined by a convergence criterion based on the L2-norm of the difference between subsequent meanflows down to a threshold 10 -8 . The closure term f is obtained by plugging u and p into the RANS equations, algebraically (see example in Fig. 2b at Re = 150).

Dataset

Data are central when considering ML based algorithms. As mentioned in the introduction, in order to obtain a model which is capable of generalization, i.e. predicting data unseen during the training process, it is necessary to collect a large amount of data. We can limit data-hungriness and overfitting problems characterizing other techniques such as the standard fully connected Neural Networks, by using a graph neural network (GNN) model detailed in §3: by exploiting the geometry and mesh independence of GNNs, we are capable of learning the underneath operator and thus overcoming the mentioned problems. We consider as input data the meanflow u and as output the stress tensor f in Eq. ( 2). The goal is to identify a data-driven model for Eq. ( 2) by means of supervised learning. We compute these quantities of interest using three meshes differently refined and running the simulations at various Reynolds numbers; in particular, the final dataset is composed of simulations ranging from Re = 40 to Re = 150 with a stride between two subsequent simulations of ∆Re = 10. Simulations at 40 ≤ Re ≤ 70 were performed on a coarse mesh (Fig. 3a), at 80 ≤ Re < 110 on a medium mesh (Fig. 3b), and at 120 ≤ Re ≤ 150 on a refined mesh (Fig. 3c). Note that we did not include Re = 110 in the training set. As shown in Fig. 3, these meshes have been refined in the most sensitive regions identified by the wavemaker (the near-wake of the cylinder) and around the obstacle boundaries, in order to increase the accuracy of the computations where the instability develops [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF]. The resulting dataset is relatively small, as it consists solely of 11 pairs of input-output flow fields. We will demonstrate the extent to which the model generalizes by performing predictions at unseen Reynolds numbers and with different geometries of the obstacle.

Methodology

The application of GNN combined with unstructured meshes revolves around the idea of graph representation of the data, which performs at its best when relations between data are relevant to the problem. In this case, we exploit the spatial distances between different points of the mesh in order to reproduce the convective/diffusive dynamics of the system and thus to learn the RANS equations operator: each node of the mesh can be coupled one-to-one to a node of the GNN, which shares with the mesh the same connection properties as its own edge features. A sketch of these connections composing a graph is shown in Fig. 4. The methodology is inherited and adapted from the work of [START_REF] Donon | Deep Statistical Solvers[END_REF], where the algorithm deep statistical solver is introduced; we refer the interested reader to this work. Here, we quickly summarize some of the fundamental features that characterize a GNN; in particular, the main process on which a GNN is based is the diffusion of the information extracted on each node and each connection across the whole network using the message passing during the updates; it can be divided in three fundamental steps:

1. Message Creation -on each node i, an embedded state associated with the vector h i is created. It is initialized with a zero state and it will embed information as long as the message passing proceeds. 2. Message Propagation -in this step, information is propagated. For unravelling the dynamics of the underneath system, i.e. producing a model for the RANS equations, we transmit the message in both directions: from each node to its neighbours and from these latter to the node itself (Fig. 4). Mathematically this is stated as φ φ φ

(k) i, j = ζ (k) i h (k-1) i , a i j , h (k-1) j , (4)
In what follows, the subscript will indicate the connections. The vector h

(k-1) i
is the embedded state from the previous layer k -1, a i j denotes the directed connections between the node i and the nodes j, while ζ (k) i is a generic differentiable operator such as multi-perceptron layers (MLP), see [START_REF] Goodfellow | Deep learning[END_REF]. 3. Message Aggregation -finally, on each node, these collected informations are aggregated, providing an updated embedded state h

(k) i for each of them h (k) i = h (k-1) i + αΨ (k) h (k-1) i , {u, Re}, φ φ φ (k) →,i , φ φ φ (k) ←,i , φ φ φ (k) ,i , (5)
where φ φ φ

(k) →,i = 1 N d ∑ N d j=1 φ φ φ (k)
j,i is the message that the node i sends to its N d neighbours; φ φ φ

(k) ←,i = 1 N d ∑ N d j=1 φ φ φ (k) i, j
are the messages that the node i receives from its N d neighbours; φ φ φ

(k) ,i = φ φ φ (k)
i,i the message that the node i sends to itself to avoid lost of information as the process advances. Ψ (k) is a generic differentiable operator that can be approximated also in this case by MLP; α is a relaxation coefficient, included in the model hyperparameters. Note that the use of the mean operator in φ φ φ

(k) →,i and φ φ φ (k)
←,i allows one to collect, for each of the i-th embedded states, the information propagating from a number of neighbours that can differ for each of the i-th nodes. Nonetheless, different permutation-invariant functions can be introduced in this step (maximum, sum, concatenation...). In our case the inputs -indicated in Fig. 4 as G i -are the meanflow and the Reynolds numbers {u, Re}, and are fed during the training step.

Following the strategy adopted in [START_REF] Donon | Deep Statistical Solvers[END_REF], we apply different trainable parameters for each layer of the GNN. At the end of the message passing process, the information on each node has been handled as described and thus carries with it information from any other point of the graph including their relative distances. In principle, the number of updates should cover the longest geodesic path that can be defined on the mesh, in order to allow the message on each node to reach any other node of the mesh; in practice, the number of updates are fixed and -in our application - optimized using genetic algorithms [START_REF] Akiba | Optuna: A next-generation hyperparameter optimization framework[END_REF]. The latest embedded state is projected back to a physical state using a decoder, once again approximated using a differentiable function that in our case is a fully connected network.

A detailed discussion of the parameters, their choice and the structure of the NNs is included in the supplementary appendix.

In short, a GNN layer is a parametrized function that takes in input a graph and returns the same graph with updated node and/or edge features, learnt during the training process and compared with the ground truth's RANS closure term f computed using DNS simulations. The resulting error leads to the definition of a loss function based on the Euclidean norm or mean square error ε, reading as

ε = k ∑ k=1 γ k-k 1 N N ∑ i=1 (x i -y i ) 2 , (6)
which needs to be minimized during the training. In Eq. ( 9), N is the number of nodes, x is a vector containing the GNN's prediction for each node, y is a vector of the same dimensions representing the ground truth, k is the number of GNN layers and γ a process parameter. Gradients of ε with respect to the weights of the GNN are computed by automatic differentiation; stochastic gradient descent is adopted for updating the weights of the net to best fit the training data (back propagation) and minimizing ε.

Results

In this section, we briefly present the results based on the GNN model discussed in §3. The implementation is based on Pytorch Geometric, a PyTorch library for building and training GNNs, while the coupling between the network and the mesh has been implemented with a Python scratch-coded script which provides an interface with the FEM library. The training of the model has been executed on 2 parallel GPU NVIDIA Tesla V100 with 32 GiB of RAM. Fig. 5 shows the training loss as a function of the number of epochs. By considering the same architecture, we checked the robustness of the training by running 10 models (Fig. 5a) on a limited number of epochs: for all cases, the drop in the loss is observed during the first 2000 epochs, as indicated by the mean and the variance. This behaviour indicates that all the trained models are characterized by comparable accuracy; we finally consider one of these 10 models and train it up to 125k epochs (Fig. 5b). Using this model, we consider in the following different test cases for assessing the quality of the prediction.

The quantitative comparison between the GNN's prediction and the DNS's ground truth from the RANS closure terms, has been performed by computing the relative error

∆f = f dns -f gnn f dns , (7)
where f dns is the closure term of the RANS equations from direct numerical simulations and f gnn is the GNN's prediction of the same term. First, we start from a proof of concept case, meaning that f has been predicted by the GNN on a case, Re = 130, that is included in the training dataset (Fig. 6a-c and Fig. 8b). In this case an error ∆f = 0.034 is As validation, we first consider cases that are not included in the dataset, by varying first the Reynolds number: in Fig. 6df and Fig. 8d the test is performed at Re = 200. This case is beyond the Reynolds number range used for training, set at Re = 150; the accuracy obtained in this case is ∆f = 0.39. In particular, we note that good agreement is achieved in the cases shown in details in Fig. 8a and Fig. 8d where the wake is fully developed, while discrepancies can be found in the immediate vicinity of the obstacle. Then, we consider changes in the configuration. In Fig. 7, the cylinder has been down shifted and still the GNN is able to show good prediction results (∆f = 0.29), confirmed in the 1D plot in Fig. 8e. Finally, we analyse a different geometry for the bluff body: we assess the prediction of the closure term on a tilted 2D square in Fig. 7. In this case a decay of the performance is observed, with an error of ∆f = 0.79. The GNN model fails at capturing the details of the areas characterized by strong velocity gradients, although we can observe an overall ability in capturing the main patterns of the spatial distribution.

Conclusions and discussions

Supervised learning is applied for the prediction of the closure of Reynolds-averaged Navier-Stokes equations. We consider graph neural networks (GNN) models, exploiting the advantages that characterize this strategy: compared to standard convolutional neural networks, the high flexibility in computing a convolution adapted on any kind of interconnected, unstructured mesh allows to extract a model from any kind of data-set. Leveraging this feature, here we consider numerical simulations based on finite elements method (FEM). This work enters a very recent line of Despite the small amount of data, we found that the model is capable of reconstructing the Reynolds stress, with a good fidelity. On one hand, these performances are possible as the model is independent from the geometry and the mesh discretization. On the other hand, this striking difference with respect to analogous approaches is due to the application of statistical learning processes inspired by the deep statistical solver algorithm discussed in [START_REF] Donon | Deep Statistical Solvers[END_REF], where the ultimate goal is the operator learning. In [START_REF] Donon | Deep Statistical Solvers[END_REF] the potentialities of this strategy were proven mathematically for unsupervised cases. Here, we demonstrate that numerically for the supervised case by concluding that -in principle -the obtained GNN based on DSS is capable of generalization.

Following this philosophy, a number of improvements and follow-ups can be introduced. First, as the modelling is strongly dependent by the quality and amount of data, while being relatively parsimonious compared to analogous techniques, methodologies that systematically account for a selection of the data can be introduced, as well as attention mechanisms weighting the most informative regions of the flow (graph attention networks). The updates of the GNN, here performed through a chain of multi perceptron layers, can be replaced by one recurrent neural network [START_REF] Hamilton | Graph representation learning[END_REF]. Finally, the terms approximated by these models could be integrated within turbulence modelling schemes and tested at higher Reynolds numbers. These research venues are currently subject to scrutiny.

A Numerical simulations

Numerical simulations were performed using a code written in FEniCS [START_REF] Alnaes | The fenics project version 1.5[END_REF]. Time resolved simulations were used for building the dataset, while the inputs and the outputs of the model were obtained by averaging these data. From the numerical viewpoint, the spatial discretization is obtained by the weak formulation based on the finite element method (FEM). In particular, the finite element used is the Taylor-Hood element, with second order elements P2 for velocity and first order elements P1 for pressure. The implemented integration scheme reads as

         3u n -4u n-1 + u n-2 2∆t + (u n-1 • ∇)u n + (u n • ∇)u n-1 -(u n-1 • ∇)u n-1 - 1 Re ∆u n + ∇p n = 0 ∇ • u n = 0, (8)
where u = (u, v) T represents the velocity vector, p the pressure and Re the Reynolds number. The time step is denoted by ∆t; the n apex indicates the values of a quantity at the current time, with n -1 the values at the previous time step and n -2 its values two time steps before. The time marching is performed by second order backward differentiation formula (BFD). The convective term is treated using the Newton and Picard methods, such that a linear system of equations is solved in each step of the temporal iteration. Boundary conditions complete the numerical setting and are reported here for completeness. Uniform flow u = (1, 0) T is imposed at the inlet. At the outlet, p = 0 is imposed, while on the far boundaries ∂ y u = 0 and v = 0. No slip conditions are considered around the obstacles; all the numerical simulations were initiated with null flow fields at t = 0. More details on the numerical scheme can be found in [START_REF] Zienkiewicz | The Finite Element Method for Fluid Dynamics[END_REF], while the validation of the code is reported in [START_REF] Guégan | Contrôle d'écouleements turbulents par apprentissage automatique[END_REF].

B Graph Neural Network details

This section focuses on the design of the graph neural network (GNN). The main code has been written in PyTorch Geometric, a library built upon PyTorch for the developing and training GNNs. Theoretical aspects and applications to RANS are described in the main article; here, we describe and comment on parameters choices and implementation.

B.1 Overall architecture

The whole framework of the GNN is described below, with reference to Fig. 9. The starting point is H 0 which represents the tensor composed by all the embedded states defined on each node, initialized at the zero state. Along with the externally injected quantities (G, the mean flow and the Reynolds number Re in our case), it is provided to the message passing process M 1 θ (see §3 of the main paper). This latter will output an updated version of the embedded state on each node, H 1 which will pass through a decoder, D 1 θ , a multi-layer perceptron (MLP) trainable function in charge of reconstructing a meaningful physical state from the embedded state (U 1 ) which, in our case, is the closure term of the RANS equations. Finally, this physical state is compared with the ground truth that comes from the DNS via a loss function. The entire process is repeated in each layer until the final layer k is reached. The last prediction (U k) represents the output of the entire algorithm. Following the intuition in [START_REF] Donon | Deep Statistical Solvers[END_REF], all the intermediate loss values are considered, in order to robustify the learning process, in a global loss function (Eq. 3.3 of the main paper)

ε = k ∑ k=1 γ k-k 1 N N ∑ i=1 (x i -y i ) 2 . ( 9
)
ε is the loss function, N is the number of nodes, x is a vector containing the GNN's prediction for each node, y is a vector of the same dimensions representing the ground truth and k is the number of layers.

B.2 Hyperparameters

Numerous parameters define the structure of the neural network, usually denoted with the term hyperparameters.

These terms are defined a-priori, before the training of the model, thus they need to be manually tuned as they can't be "learnt" during the training process. In the following we list some of them, by making a distinction between model hyperparameters ( §B.2.1) and process hyperparameters ( §B.2.2).

B.2.1 Model hyperparameters

A model hyperparameter defines the capacity of the neural network, i.e. the ability of the model to represent functions of high complexity. Thus, the capacity is directly related to the possibility of approximating a large variety of nonlinear functions. Two parameters directly impact on the capacity of the considered GNN, the embedded dimension and the number of updates.

Embedded dimension -this hyperparameter defines the dimension of the embedded vector encoding the observed quantities of the physical system. During the training process, these encoded observables are injected at each layer, whose number is defined by the number of updates k. Note that the embedded state does not have a direct physical meaning.

Number of updates, k -the number of updates of each embedded state is denoted by the variable k; the GNN behaviour relies on the message passing process (see §3 in the main document): embedded states assigned to each node are updated at each step of the message passing iteration in order to acquire information from each of the neighbouring nodes while the process progresses.

B.2.2 Process hyperparameters

The training is defined by the process hyperparameters. Tuning these hyperparameters deeply modifies the duration of the training, its computational costs and the way the weights are adjusted while the model evolves.

Update weight, α -the embedded state's update weight α is a relaxation parameter that allows to scale each update of the embedded states with the previous one during the message passing loop, as stated in Eq. (3.2) of the paper and reported here for clarity of exposition

h (k) i = h (k-1) i + αΨ (k) h (k-1) i , {u, Re}, φ φ φ (k) →,i , φ φ φ (k) ←,i , φ φ φ (k) ,i . (10)
Here, for each node, h

i is the embedded state at the current layer k; h (k-1) i the one at the previous layer; Ψ (k) is a generic differentiable operator such as MLP; G i is the externally injected input. Moreover, we recall that φ φ φ (k) →,i is the message that the node i sends to its neighbours, φ φ φ (k) ←,i are the messages that the node i receives from its neighbours, φ φ φ (k) ,i is the message that the node i sends to itself to avoid losing information as the process moves on in time.

Loss function weight, γ -this hyperparameter is used to control the importance of the contributions associated with each update of the embedded states, as reported in the Eq. (3.3) of the paper and in Eq. (2.1) of the present appendix. The relative importance of the message is assigned through an exponential increasing function, thus meaning that the latest steps, which are supposed to be the richest in information, will have the highest importance in the process.

Optimizer -the training algorithm consists of an optimization, where the objective is to minimize the cost function in Eq. (3.3). One of the most used in scientific machine learning is the Adam Optimizer, used in the present work and based on the stochastic gradient descent (SGD) method Diederik P. Kingma [2017]. The length of the step taken by the SGD method in the steepest descent direction is defined by the learning rate.

Learning rate, LR -the optimization underneath the learning process requires the definition of a LR; to this end, the scheduler ReduceLROnPlateau has been used to let this parameter to be dynamically adjusted, in order to control the SGD method online during the training process. This hyperparameter alleviates the risk of having a learning rate being too large in vicinity of the convex point of local minima, or too small, thus leading to, respectively divergence or slow convergence.

B.2.3 Optimizing the hyperparameters

The GNN architecture presented in the paper is defined upon a certain number of hyperparameters. This does not pose problems in deep learning when the expressivity of the NN (e.g. number of neurons, number of layers, etc.) is enough to represent the complexity of the problem at play. However, since the GNN is trained to fulfill a specific task and because the computational cost of the GNN inference has to be affordable compared to the most sophisticated turbulence models available today, we tried to keep the GNN as parsimonious as possible. To that end, a further optimization is required on the set of hyperparameters defining the architecture. Unfortunately, standard gradient based optimizers cannot be employed when dealing with integer numbers (i.e. number of neurons, number of layers). For this purpose gradient-free algorithms have been used. There exists dedicated libraries that can automate the tuning process through all the possible sets of hyperparameters by trying and appropriately pruning the unpromising sets of them. In this work we apply the library Optuna [START_REF] Akiba | Optuna: A next-generation hyperparameter optimization framework[END_REF], an open-source package which combines efficient searching and pruning algorithms. By exploring the complex solution space, a number of optimized combinations of hyperparameters is found, and the one that outperforms the others on the monitored validation metrics is characterized by the following set 1. Embedded dimension, 18 2. Number of GNN layers, k = 87 3. Update relaxation weight, α = 0.01 4. Loss function weight, γ = 0.1 5. Learning rate, LR = 3 × 10 -3 , as maximal/starting value.

B.3 Boundary conditions of the GNN

A specific treatment for the nodes corresponding to the boundary conditions has been implemented. We start by recalling that during the message passing iterations, each node provides and retrieves informations contained from the neighbouring nodes. This process stems from the bi-directionality of the connections a i j between each node. On the nodes corresponding to Dirichlet boundary conditions, namely the no-slip condition around the obstacle, we can apply a different strategy: in order to keep fixed the value on the boundaries, the direction of message propagation is only outward. In this way the BCs nodes are still able to provide information to their neighboring nodes but, since they can't receive anything back, they will be trained by keeping values equal to the enforced initial boundary conditions.

Figure 1 :

 1 Figure1: Sketch of the computational domain. The cylinder with unitary diameter is considered as a reference case. Boundary conditions can be found in the text.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: (a) Streamwise component of the meanflow and vorticity isolines for the flow past a cylinder at Re = 150. (b) For the same case, the streamwise closure term of RANS equations Eq. (2) is shown. In both cases, only a portion of the domain is shown.

Figure 5 :

 5 Figure 5: Training curve of 10 different models: mean loss (blue line), together with the standard deviation (-σ , +σ ) (coloured zone), (a). The green line indicate the actual model used in this work. The latter is shown for the final training performed over 125k epochs in (b) as a green line and compared to the validation loss (yellow line).

Figure 6 :

 6 Figure 6: Cylinder case. Proof of concept at Re = 130: closure term from DNS (a), GNN prediction (b), difference (c). Unseen case, Re = 200: DNS (d), prediction (e), difference ( f ).

Figure 7 :

 7 Figure 7: Unseen cases, downshifted at Re = 110: closure term from DNS (a), GNN prediction (b). Tilted square at Re = 110: DNS (c), prediction (d).

  shown in Fig.6c. In Fig.8b, the same case is shown as a function of y in 7 different locations along x behind the cylinder in the plot.

Figure 8 :

 8 Figure 8: Reynolds stress profiles along y at 7 different locations past a cylinder: comparison between ground truth (red solid line) and GNN predictions (gray dotted line). From left to right: (a) Re = 110, unseen case; (b) Re = 130; (c) Re = 150; (d) Re = 180, unseen case; (e) downshifted cylinder at Re = 110, unseen case.

Figure 9 :

 9 Figure 9: Fig. inherited from the paper of Donon et al. [2020] showing the overall framework of the GNN.

Acknowledgements. The authors acknowledge M. Nastorg (LISN) for discussions. Funding. The Ph.D. fellowship of M. Quattromini is supported by the Italian Ministry of University. A part of the research was funded by the grant PRIN2017-LUBRI-SMOOTH of the Italian Ministry of Research and ANR-21-REASON from the French Agency for National Research.


Malika Boualleg
email: malika.boualleg@ifpen.fr

Dr Dorothée Laurenti
email: dorothee.laurenti@ircelyon.univ-lyon1.fr

Dr Nathalie Tanchoux
email: nathalie.tanchoux@enscm.fr

FCCat 3 in Ronce-les-Bains: First Post-Covid catalytic gathering for the French Catalytic Community

we report here on the organization and proceedings of FCCat 3, the 3rd French Conference on Catalysis, and the first in person since GECat 2019 (the annual conference of the French catalysis community) and the covid pandemic, which took place in Ronce-les-Bains from 31 May to 3 June 2022. This international conference is held every 3 years since its first edition in 2016 and aims to bring together all catalysis communities, geographically or thematically.

Conference report

After two years of Covid pandemic issues and its consequences, which were the cancelation of the GECAT 2020 and a virtual edition of GECAT 2021, we were glad to organize the third edition of the French Conference on Catalysis (FCCAT) which took place at Ronce les Bains, Charente Maritime, from May 31 th to June 3 rd 2022. It was also a great success and the joyful occasion to gather people who haven't meet each other for two years. During these past two years, people have continued collaborations, meetings and conferences, but only…virtually and it was not easy to meet new people to establish collaborations or even just share a nice moment together. This organization of this conference continues the long history of conferences of the very dynamic catalysis community in France.

Since decades, numerous French scientists from industrial and academic laboratories played a major role on Catalysis and in the development of Catalysis at national and international levels. To promote the communication and the collaborations between French scientists, the GeCat, (Groupe d'Etude en Catalyse -Study group in Catalysis), was created in 1970, and has been since then in charge of organizing the national conference gathering every year the whole French Catalysis Community. After years of being a thematic group belonging to the SCF, The GECat has been now fully linked to the Catalysis Division of French Chemical Society (Groupe d'Etude en Catalyse -Study group in Catalysis) since 2020.

The French Catalysis Division (DivCat) was founded in 1983 and is the one of the scientific entities of the French Chemical Society (SCF). Its purpose is to promote excellence in the field of Catalysis at national, european (via the European Federation of Catalysis Societies, EFCATS) and international levels (via the International Association of Catalysis Societies, IACS). The DivCat also favors the exchanges between academia and industry in all fields of catalysis, and the opening toward the others scientific fields.

French conferences are thus annually organized by the French catalysis community to exchange ideas and scientifics insights. Among its, the most recurrent are the GECAT and the FCCat.

The meeting comprises invited plenary and keynote lectures and selected oral and posters communications. In 2014, the DivCat together with the GECat executive board proposed to make a new scientific event, more ambitious, in the form of an international conference organized by the French catalysis community: the French Conference on Catalysis (FCCat). The Groupe Français des zéolithes (GFZ, French Groupe on Zeolites), as well as the Industrial Chemistry Division, the Coordination Chemistry Division and the Organic Chemistry Division of SCF were also part of the organization. The very first edition FCCat took place in Fréjus, French Riviera, in May 2016. FCCat1 addressed a large spectrum of topics, from homogeneous to heterogeneous catalysis, from fundamental to applied questions. About 240 participants joined, which appeared to be the ideal scale for very intense scientific discussions in a convivial atmosphere. 35 % of non-French participants were present, and the industrial participation was very strong (22 %), promoted by the "industry half-day" devoted to relationships between industry and young researcher.

After the success of the first edition, it has been decided to continue to organize this scientific event, on a 3-year periodicity, the FCCat replacing the GECat in these years. The second FCCAT edition took also place in Fréjus, in June 2019 and gathered more than 210 attendees. 30 % of non-French participants were present, and the industrial participation was again very strong (20 %), promoted by another "industry half-day" devoted to relationships between industry and young researchers. As usual and as quickly became the tradition, heterogeneous, homogeneous, hybrid catalysis researches were presented at the second FCCAT edition.

The third edition took place this year in May-June 2022 and gathered 180 attendees in a wonderful place by the sea in Charente Maritime. The number of attendees was great, especially considering that the pandemic still prevented all foreign scientists to join. The FCCat3 conference started with an international panel of speakers : Professor Susannah Scott coming from Santa Barbara University (US) presented the catalytic valorization of renewables or recycled feedstocks for circular carbon economy. Her presentation was focused on bio-sourced molecules and plastic catalytic transformations with particular attention paid to the effects of solvents in these reactions. Just after that plenary conference, Professor Jeremy Luterbacher (keynote speaker) from EPFL (Lausanne, Switzerland) spoke about the simplification and the control of catalytic upgrading of lignocellulosic derivatives, thanks to chemical functionalization to prepare bio-sourced materials. The conference was marked from the very beginning by the need to use catalysis to improve our carbon footprint and valorize bio-resources!

In order to achieve the catalytic conversion of bio-resources and small molecules, new catalytic performant systems are always required. In this spirit, the next day, Professor Damien Debecker (Univ. Louvain, Belgium) further presented, in a very inspiring plenary lecture, the preparation of supported nanocatalysts with innovative techniques with applications for CO 2 and glycerol transformation.

The second day, oral flash presentations were presented on various topics, and, then, we divided in 2 parallel sessions: one parallel session was dedicated to porous catalytic materials, while the second one was focused on catalytic sugars and cellulose conversion.

This 31 st of may morning finished with a keynote of Professor Stéphane Baudron (Univ. Strasbourg) who presented deep eutectic solvents effect on MOF preparation and materials properties. In the same afternoon, we had again parallel sessions with one devoted on atom economy in synthetic chemistry and one about H 2 production. Professor Moniek Tromp (Groningen Univ., The Netherlands), who had to cancel her venue due to flight cancellation but was keen to present by visio conference her keynote on insights in dynamic structures and electronics for catalysis. Then, before the poster session that led us to diner, we had a short lecture of Charlotte Gers Panther from Chemistry Europe, who invited all of us to submit a communication to this ChemCatChem special collection.

The second day was also rich, beginning with a plenary lecture from Doctor Thibault Cantat (LCMCE, CEA Univ. Paris Saclay) on the catalytic strategies for the reduction of C-O bonds with applications towards CO 2 and other small molecules, like formic acid or CO, valorization. Then the full assistance attended the second flash session of the conference before dividing into the two parallel sessions on CO 2 valorization and metal sulfide catalysts/(hydro)desulfurization. The morning ended with the keynote of Professor Juan Ramon Gonzalez Velasco (Univ. Basque Country, Spain) who spoke about multifunctional materials for CO 2 capture and conversion towards synthetic natural gas. The afternoon was dedicated to social activities: we privatized the "train des mouettes" which included a guided pedestrian visit of Mornac sur Seudre, labelled "most beautiful village of France" and the guided visit of the salt marshes. The second poster session took place in the evening with a wonderful tasting of Oysters of Mornac before diner. Thursday, 2 nd of June, began by a very didactic and stimulating plenary talk by Doctor Florence Epron (IC2MP University of Poitiers). She showed all her expertise on bimetallic catalysts synthesis, especially by innovative methods developed in her group, like redox reduction-oxidation reactions, which allow the deposition of a second metal at the surface of monometallic supported nanoparticles. This method leads to a better control of metals deposition and particle size. These catalysts have wide fields of application in gas and liquid-phase catalysis.

As we wanted to mix all the communities interested in catalysis and ensure a wide audience, we then We then organized a very successful industrial session and round table destined to young (and less young!) people, who might have to look for an employment in the near future. We are very grateful to Benjamin Faure (Solvay), Patrick N'Guyen (Saint-Gobain), Fabien Grasset (MANE), Xavier Baucherel (Jonhson Matthey), Stefan Schunk (BASF/HTE), Hélori Salambier (Arkema) and Juliana Aparicio (Dielix/Sarpi Veolia), who came to participate to this industrial session. They all exposed the specific activity of their company and shared their views on the job market and opportunities in their company. Students from the assembly could then ask questions which led to a lively and enriching debate. We hope that it gave many students good tips and ideas for the continuation of their career.

The two parallel sessions of this day were devoted to metals and carbides catalysts on one side and to Hydrodeoxygenation and CO 2 and CH 4 transformation on the other side, and both sessions led to stimulating discussions and insights into the research projects carried out in French or European laboratories.

As is traditional at the French Catalysis conference every year (GECat or FCCat), we invited the three After a nice gala dinner and musical evening, the last plenary lecture was given by Guylène Costentin, who showed us very nice examples of the synthesis and use of hydroxyapatites for heterogeneous catalysis. The last two parallel sessions of the conference were devoted to bifuntional and acid-base catalysis, as well as catalysis for depollution.

During the conference, several juries awarded oral and poster prizes to PhD and post-docs. Two best oral presentations, sponsored by IFPEN and NJC, were awarded, 2 best flash presentations (sponsored by NJC and IFPEN) were also awarded. And finally, 3 best posters presentation awards were given by Chemistry Europe.

During the whole conference, we had a wonderful experience with the "NanoRêVe" device thanks to Pascal Roland and Professor Arnaud Travert from LCS (Caen). "NanoRêVe" is an immersive virtual reality platform that allows people to explore objects on a molecular scale and proposed us an immersive experience with some of our favorite catalysts (https://nanoreve.cnrs.fr/)! Many instrument companies (Vinci Technologies, Anton Paar, Thermoscientific, Equilabo, Ablatom, Top Industrie, Magritek, Micromeritics and CAS) were also present to show us their new products and discuss with participants during all coffee breaks and more. We warmly thank them for their presence. Finally, our thanks go to the whole organizing committee who worked hard for this conference to happen, and to all the participants who came to present results, discuss science, eat, drink, socialize, and made this conference a success.

Figure 1 .

 1 Figure 1. "Train des mouettes" used by the participants during the recreational afternoon

  programmed a keynote lecture by Doctor Alessandro Fortunelli (CNR, University of Parma), a theoretical chemist, who gave us an exciting lecture on combined theoretical/experimental approaches to understand the parameters affecting the performances of Pt-based catalytic systems for Oxygen Reduction Reactions. It is a nice example of the interest of mixing theoretical and experimental data to understand how a catalyst work and help design new catalysts.

  DivCat prizes winner to give a lecture about their work. The first lecturer, Gabriel Durin, winner of the best PhD Award, sponsored by BASF, gave an enlightening talk about his PhD project supervised by Thibault Cantat (CEA) on the catalytic hydrogenolysis of halosilanes to give hydrosilanes. The Young Researcher award, sponsored by IFPEN, was awarded to Prince Amaniampong, a young CNRS researcher, newly recruited at IC2MP, who explained his nice projects on catalysis under ultrasound irradiations, mainly applied to the transformation of biomass. The third prize, which was newly created this year by the DivCat, is the Innovation award, and it was awarded to Bertrand Guichard, Engineer and researcher at IFPEN. He talked about the development over the years of new classes of catalysts for hydroprocessing and hydrocracking.For the past few years, all the French catalysis community has been also involved in the preparation of the next International Conference on Catalysis (ICC Lyon 2024). Céline Chizallet (IFPEN) presented us the latest steps of the organization to keep the community informed and motivate everyone's implication to make this event dynamic. Effectively, for the first time, the worldwide catalysis community will have a participatory and interactive opportunity to contribute with its ideas and energy to the preparation of an ICC. So have a look to build your ICC and ICC2024 web sites (https://www.icc-lyon2024.fr/ & https://build-your-icc.org/ )

Figure 2 .

 2 Figure 2. A participant testing the "NanorêVe" VR platform by walking inside a zeolite pore

Figure 3 .

 3 Figure 3. Picture of all the happy participants!

  Indeed, every year since 1970, the GECat executive board composed by 10 members of the main French laboratories in, mostly heterogeneous, catalysis, organizes a conference gathering the

	French-speaking	catalysis	community	around	three	different	specific	themes
	(https://new.societechimiquedefrance.fr/divisions/catalyse/gecat/), where young researchers'
	contributions are highly encouraged and favored. A different and very convivial venue is chosen each
	year all around France which gathers around 100 participants.

Table of content entry Report

on FCCat 3, French Conference on Catalysis, which took place from May 31 st to June 3 rd 2022 in Ronce-les-Bains, France Authors Biography

Dr. Dorothée Laurenti Dorothée Laurenti obtained a PhD from University of Aix-Marseille in 1999 under the guidance of Pr M. Santelli in Organic Synthesis and homogeneous catalysis. Then she moved to York University (UK) and work in collaboration with Pr J.H. Clark (Green Chemistry Center) and Rhodia company on industrial process with supported homogeneous catalysts. Afterwards, she developed more expertise in catalysis focusing on plastic transformation with catalysts prepared by surface organometallic chemistry under the supervision of Dr J.-M. Basset and E. Kuntz in Lyon. In 2002, she joined the Institute of Research on Catalysis and Environment in Lyon (IRCELYON) as CNRS Researcher to work with M. Vrinat and C. Geantet on heterogeneous catalysis for hydrotreatment of fossil fuels. Her main research interest focuses now on biomass and wastes valorization by heterogeneous catalysis towards biofuels and bio-chemicals. She leads the group CATREN at IRCELYON and is the director of the French group of research on Thermochemical Conversion of Biomass and Wastes (GDR Thermobio). She served as GECAT board member (2010-2012), DivCat board as secretary (from 2019) and is one of the executive organizers of the next ICC2024 conference that will take place in Lyon next July.

Web site: https://www.ircelyon.univ-lyon1.fr/syrcel/card/DLA/


Anthony Bernard
email: anthony.bernard@inra.fr

Fabrice Lheureux

Elisabeth Dirlewanger

SSR genetic diversity assessment of French walnut

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

➢

  Genetic diversity and population structure assessment of INRA's walnut (Juglans spp.) germplasm collection using SSR markers Study of Juglans spp. genetic resources for the implementation of a marker-assisted selection Context • Important economic development of walnuts • France: insufficient choice of cultivars • Global competition, climate change, European plan Ecophyto II • Genetic improvement: important lever SSR Genetic Diversity Assessment of the INRA's Walnut Germplasm Collection INNOV'noyer project INNOV'noyer project (2) ➢ Study of Juglans spp. genetic resources for the implementation of a marker-assisted selection Goals • Initiation of a new breeding program • Basic research: -Genetic diversity evaluation of germplasm repository -Genetic determinism identification of traits of interest • Applied research: establishment of necessary tools for marker-assisted selection achievement SSR Genetic Diversity Assessment of the INRA's Walnut Germplasm Collection Action plan • Characterization of the 253 accessions from Genetic Resources Center Juglans: -Available phenotyping data (from 1970 to 2006) and new phenotyping in 2017, 2018, 2019 -Traits: tree, phenology, fruit and diseases • Genotyping: SSRs and 600K SNP array Axiom® • Genome-Wide Association Study 3 INNOV'noyer project INNOV'noyer project (3) ➢ Study of Juglans spp. genetic resources for the implementation of a marker-assisted selection

236 236 236 236 236 236 236 236 236 236 236 236 236 236 236 236 236

		SSR Genetic Diversity Assessment of the INRA's SSR Genetic Diversity Assessment of the INRA's SSR Genetic Diversity Assessment of the INRA's
		Walnut Germplasm Collection Walnut Germplasm Collection Walnut Germplasm Collection Walnut Germplasm Collection
	36 of related species GRC Juglans -Fruit Tree Experimental Unit (INRA) 253 accessions Genetic diversity 253 accessions + 15 SSR markers WGA 376 Genetic diversity of J. regia collection Bar plot K=2 600K SNP genotyping WGA 001 Locus A H o H e F IS WGA 001 0.65 0.71 0.08 on 170 accessions WGA 001 WGA 004 WGA 009 WGA 027 WGA 069 WGA 072 WGA 202 WGA 276 WGA 349 WGA 376 JR 0160 JR 1739 JR 1817 JR 6160 JR 6439 Genomic SSRs from J. nigra -Dangl et al., 2005 -Woeste et al., 2002 EST-SSR from J. regia -Dang et al., 2016 WGA 004 0.46 0.50 0.08 WGA 009 0.59 0.66 0.11 H e > H o : deficiency of sieboldia na 234 234 234 236 236 238 238 238 238 238 240 240 242 242 242 242 242 242 242 242 244 244 244 244 244 244 246 248 248 248 250 250 250 255 255 255 255 257 257 257 259 259 259 261 261 261 263 263 263 265 265 265 269 269 269 271 271 271 273 sieboldia na 184 184 184 195 197 198 198 198 199 199 199 199 200 200 200 201 201 201 201 201 202 202 202 203 203 203 203 203 203 203 203 203 204 205 205 205 207 207 207 207 208 208 208 213 213 213 219 219 Max 0.68 0.78 0.37 Total 609,658 100 Off-Target Variant Call Rate Below Threshold He, Ho, FIS lower for Population 2 210 210 210 Min 0.10 0.15 Mean 8.92 0.47 0.56 0.17 Pairwise FST = 0.101 Other 66,941 10.98 0.07 OffTargetVariant 4,787 0.79 206 206 206 WGA 027 0.43 0.46 0.07 WGA 069 0.57 0.76 0.25 WGA 072 0.32 0.42 0.24 WGA 202 0.68 0.78 0.13 WGA 276 0.63 0.76 heterozygotes JR 1817 and JR 6439: lower Category Number of Poly High Resolution No Minor Homozygote Mono High Resolution % of markers markers 0.17 WGA 349 0.49 0.78 0.37 WGA 376 0.61 0.66 H o and H e PolyHighResolution 397,921 65.27 0.08 JR 1817 0.10 0.15 0.33 JR 6160 0.43 0.49 0.12 JR 6439 0.14 0.17 0.18 NoMinorHom 75,564 12.39 Population 1: Eastern European and Asian accessions (x 63) MonoHighResolution 36,684 6.02 473,485 SNP -Population 2: Western European and American accessions (x 127) CallRateBelowThreshold 27,761 4.55 77.66%
	5 10	Genetic diversity and population structure (1) Genetic diversity and population structure (4)

SSR Genetic Diversity Assessment of the INRA's Walnut Germplasm Collection INRA's walnut germplasm INRA's walnut germplasm 217 of J. regia L. Genetic diversity and population structure (2) SSR Genetic Diversity Assessment of the INRA's Walnut Germplasm Collection Genetic diversity and population structure (3) SSR Genetic Diversity Assessment of the INRA's Walnut Germplasm Collection Ancestry coefficient Genetic diversity and population structure (6) SSR Genetic Diversity Assessment of the INRA's Genome-wide SNP genotyping (1) SSR Genetic Diversity Assessment of the INRA's Walnut Germplasm Collection Genome-wide SNP genotyping (2)

Acknowledgements Acknowledgements ➢ Eric Germain, former head of INRA's breeding program from 1977 to 2007 ➢ Tree Experimental Unit in Toulenne and Prunus/Juglans Genetic Resources Center ➢ « Région Nouvelle-Aquitaine » for funding the project ➢ « ANRT » for funding the PhD ➢ « BioGEVES » laboratory for SSR amplifications

Step 1 20 runs: K from 1 to 10 Burn-in period: 5,000 Markov Chain Monte Carlo: 50,000 Admixture model Correlated allele frequencies Plateau criterion + ΔK method K=2 ?

Step 2 10 runs: K from 1 to 5 Burn-in period: 100,000 Markov Chain Monte Carlo: 750,000