Solving the Uniform Density Constraint in a Stochastic Downscaling Model

The present work aims to contribute to the development of a numerical method to compute small scale phenomena in atmospheric models, getting rid of any mesh refinement. In a domain, typically a mesh of a numerical weather prediction model, we simulate some particles that are moved thanks to a Stochastic Lagrangian model adapted from the PDF methods proposed by S.B. Pope. We estimate the Eulerian values of the required fields, thanks to the computation of a local mean value over an ensemble of particles. We are thus using a stochastic particle method. At small scale, our atmospheric model imposes that the mass density is constant in the domain. As a consequence, the particles have to be uniformly distributed at every time step of the particle method. We aim to use D.P. Bertsekas Auction algorithm in order to satisfy this constraint. Naturally, the transport cost will have to be minimum. This is a problem of 3D optimal transport, which is known to be difficult.


Introduction
This work aims to contribute to the development of a numerical tool, called SDM (Stochastic Downscaling Method), that was initially proposed in [19]. It consists in introducing a new method for the simulation of the wind at small scales. The large scale behavior of the ground wind is known to be of major importance for our society, in particular in the context of global warming. But the local wind variability is also a key mechanism, for instance when one wishes to build a bridge, or optimize the location of a wind turbine. Ground winds also vary at small time and space scales, that are unfortunately far smaller than those solved by classical weather/climate forecasting models. Thus, there is a need to imagine and propose some new numerical techniques that would ensure a satisfying prediction of the wind at sufficiently small scales: this is the general objective of the downscaling techniques. Whereas deterministic refinement methods, based on a posteriori estimates [20], are already well developed, some of the authors introduced SDM [19], an algorithm of a totally new type in the framework of downscaling methods, based on a stochastic Lagrangian model, solved by a particle method (see Section 1.1 below). Indeed, given the coarse prediction provided by a Numerical Weather Prediction model, SDM simulates a set of fluid particles in a chosen domain of computation, in order to refine the numerical prediction over it. For the time being, the vertical size of the SDM domain is sufficiently small to consider a constant mass density. This implies that the numerical particles must satisfy the uniform density constraint in the domain. At each time step, the rearrangement of the particles under minimal displacement is the most expensive procedure, in term of number of operations. To this aim, we evaluate in this paper the performances of the Auction Algorithm, introduced in [6], and its ability to be coupled with SDM. At this time, SDM uses methods expected to be nearly optimal whereas the Auction Algorithm is optimal. This tool will be our reference algorithm; in this paper, we validate the heuristic parameters of Auction Algorithm, and test its complexity on several configurations of particles.
This article is organized as follows: we start in Section 1.1 with a general description of Lagrangian stochastic models, before introducing in 1.2 the constraint of uniform density, required by the model. Then, in Section 2, we present the Auction Algorithm. Finally, Section 3 provides some numerical results, and analyzes the claimed properties of the Auction Algorithm thanks to different experiments.

Stochastic Lagrangian Models
We now give a general description of Lagrangian stochastic models, in order to introduce the specific problem tackled in this paper, that is the uniform density constraint. For more details, the reader is referred to [16,17] and the related bibliography.
We are interested in the behavior of a turbulent incompressible fluid in an open set D of R 3 , with constant mass density ρ. In the statistical description of turbulence, all the physical fields (here the velocity U and the pressure P) are considered as random fields: they can be written, following the so-called Reynolds decomposition, as the sum of a deterministic part and a random turbulent part, namely: Applying the Reynolds operator to the incompressible Navier Stokes equations, we obtain the so-called Reynolds Averaged Navier Stokes (RANS) equations: where U (i) stands for the ith component of U , with i = 1..3, and P is solution of the following Poisson equation: The matrix { u (i) u (j) } 1≤i,j≤3 , called the Reynolds tensor, is defined by: The system (1.1) is not closed, and hence one needs to provide additional information on the Reynolds tensor, in order to close the system and compute some solution. Different ways of closing system (1.1) have been introduced over the years: for example k − ε closure, or second order Rotta closure, etc. The interested reader is referred to [14] for further details on the closure of the RANS equations. If we denote by f E (t, x; ·) the law of the random variable U (t, x, ω), it is to be noticed that a closure on f E is equivalent to the RANS closure.
From 1985, S.B. Pope started to introduce some completely different way of closing the RANS equations. His idea consists in describing the fluid as an ensemble of Lagrangian particles, introducing a stochastic model, whose state variable is (X t , U t , φ t ), that respectively correspond to the position, the velocity and, possibly, some other physical quantities attached to a fluid particle, such as a temperature or a turbulent frequency. The dynamics of the fluid particle is ruled by a system of stochastic differential equations (SDEs) of the form where (W t , t ∈ R + ) and (W t , t ∈ R + ) are two independent Brownian motions. The system (1.3) intends to model f E and the laws governing the Eulerian fields U (t, x, ω) and φ(t, x, ω) in the sense that U t = U (t, X t , ω) and φ t = φ(t, X t , ω).
The design of the coefficients D * and B * (the subscript * holds for U and φ) in equations (1.3b) and (1.3c) determines the modelling context and its closure: for SDM, their expressions are given in [3]; for engineering applications such as turbulent fluids in combustion, or models including the temperature, see [10,11]. These coefficients D * and B * depend on statistics of the state variables (X t , U t , φ t ) (for instance the mean velocity, the turbulent kinetic energy, the dissipation), according to the model closure. Hence those coefficients D * , B * are non linear in the sense of McKean [12].
Assume the existence of the density f L (t; x, V, θ) for the Lagrangian variables, then at each time t, the measure f L (t; x, V, θ)dxdV dθ is the law of (X t , U t , φ t ) solution of (1.3). Considering an incompressible fluid, with constant density ρ, then the link between f E (t, x; V, θ) and f L (t; x, V, θ) (see e.g. [15] for details) is: which allows to compute any mean fields required in D * and B * by the formula: Moreover, and to end this part, a fluid particle satisfying (1.3) and (1.2) also satisfies [8] : The pressure term in (1.3b) produces these two equations (1.5). In SDM, we simulate the equations (1.3) without the pressure term (prediction step), and then correct the fluid density and the mean velocity according to (1.5). Since the coefficients D * and B * depend on the law governing the fluid particle in a non-local way (mean field interaction), a stochastic particle method is required to discretize (1.3). Next section firstly introduces the Particle in Cell method, and the splitting scheme (prediction + correction). More particularly, the aim of this work is to present a method to study and solve the constraint (1.5a), called the uniform density constraint.

Uniform Density Constraint
In a parallelepipedic domain D, we present the discretization scheme to solve the splitted problem: (1) The dynamics of each particle k is characterized at time t n = t 0 + n∆t by the discrete variables (X n k , U n k , φ n k ). The coefficients D * and B * are non-linear, and depend on each particle k, and on the statistics given by the N particles.
(2) Prediction step: solve the discretized version of (1.3) without the pressure term (3) Correction step: modify the state variables in order to satisfy (1.5).
In the context of building the numerical tool SDM that was mentioned above, we have to deal with a parallelepipedic domain D.
In this domain, we describe equations (1.3) (without the pressure term) thanks to a system of N interacting particles. The dynamics of each particle k is characterized at time t n = t 0 + n∆t by the discrete variables (X n k , U n k , φ n k ). The operators D U and B U of equation (1.3b) are then determined at time t n by D n ; the same holds for operators D n φ and B n φ of equation (1.3c). The discretized system at time t n+1 is the following : } using a local approximation, as in the Particle in Cell (PIC) method, see [18]. Among the available PIC interpolations, we use the so-called Nearest Grid Point (NGP) method, which consists in dividing the domain D into cells (see Figure 1), and in evaluating the Eulerian quantities, as well as the fluid density, at a point x ∈ D by a Monte-Carlo approximation which uses all the particles k belonging to the cell C, x ∈ C. After proceeding to the advancement step (1.6), the equations (1.5) have to be fulfilled by the system of particles. We focus here on the the uniform density constraint (1.5a), which consists in ensuring that each cell of the discretized domain (see Figure 1) has the same number of particles. To ensure a good convergence of the approximation in each cell, the number of particles in each cell N pc has to be important (greater than 1000). That is, given a set of N particles of positions (X i ) 1≤i≤N , we have to move them to some new positions (Q j ) 1≤j≤N , such that every cell of the domain contains N pc particles. The choice of final positions (Q j ) 1≤j≤N is obviously not unique, but they have to follow the uniform law: in SDM, the positions (Q j ) 1≤j≤N are constructed using the uniform law in each cell, such that all the cells contain exactly N pc particles.
In addition, we use the physical properties of the Lagrangian particles (U k , φ k ) to compute the Eulerian fields (e.g. U n ) (see [19] and Equation (1.4) above). Consequently, the particles should not move too far from their initial position during the transportation problem, in order to maintain locally the physical information they transport.
That is why we need the transport cost from (X i ) 1≤i≤N to (Q j ) 1≤j≤N to be minimum; if we define the distance between an initial position and a final one: where | · | is the Euclidean norm of R 3 , then we are faced with a classical problem of optimal transport: where S N stands for the set of permutations over {1, . . . , N }. More explicitly, we want to find a feasible assignment, i.e., a set of pairs {(X 1 , Q j(1) ), . . . , (X N , Q j(N ) )} such that j(1), . . . , j(N ) are all distinct, that minimizes the overall cost . This is an instance of the classical assignment problem appearing in network optimization. Both from the continuous and discrete viewpoints, the optimal transport problem in dimension 3 is known to be challenging (see [2,9,13]), whereas in dimension 1, it simply consists (in the discrete case) to sort the particles [2], with a computational cost of O(N log(N )), where N is the number of particles. The assignment problem (P ) can be understood as a search, among all the N ! permutations, of one permutation (possibly not unique) satisfying (1.7). Obviously, the aim is to solve the problem using a method with an efficient computational cost rather than factorial. In this direction, several algorithms have been developed and the general upper bound O(N 3 ) for the complexity was shown (see the simplex method of [1] or the so-called Hungarian method). The method we use to reach the solution of our particular assignment problem is the Auction Algorithm developed by D.P. Bertsekas, since its practical performance is competitive with stateof-the-art codes (see [6] for a short presentation). The aim of this work is to study the Auction Algorithm in our case and its complexity (we call complexity the total number of unitary assignments (X i , Q j ) during the auction process), and to see in which conditions it is an interesting algorithm for our application.

Correction of Particle Positions
In this section, we present the Auction Algorithm and some of its improvements, introduced by Bertsekas.

Naive Auction Algorithm
Having the form of a maximization problem, the classical assignment problem deals with benefits rather than costs. More precisely, there are N persons and N objects that we have to match on a one-to-one basis supposing that there is a benefit a ij for matching person i with object j. The aim is to assign persons to objects so as to maximize the total benefit (see [6]). Mathematically, the problem is: (2.1) Setting by d ij = −a ij , i.e., passing from benefits to costs, the assignment problem (2.1) is equivalent to our problem (1.7). Suppose that position Q j has a price p j in the sense that we must pay p j so as one particle can acquire this position. Thus, we define the total cost for moving one particle from initial position X i to final position Q j as where d ij is the transport cost. The algorithm seeks to assign each X i to some Q j(i) with minimum total cost, i.e., for every i ∈ {1, . . . , N } find j(i) ∈ {1, . . . , N } 3) The algorithm starts with the empty assignment and the zero price vector. At the beginning of each iteration there are a partial assignment and a price vector. The algorithm is applied until all (X i ) 1≤i≤N become assigned.
Step 1. Select randomly an unassigned initial position X i .
Step 3. Compute the bidding increment Step 4. Increase the price of position Q j(i) by γ i and assign Q j(i) to X i .
Step 5. The initial position that were assigned to Q j(i) at Step 2 becomes unassigned and then go to Step 1.
Each iteration of the algorithm can be seen as a real auction where the bidder X i searches for his preferred object Q j(i) and the second best object. Then X i raises the price of Q j(i) by γ i since at this stage Q j(i) rests his preferred object but it becomes less attractive to the others bidders.

Auction Process
Unfortunately, as shown in [5,6], the Naive Auction Algorithm (Algorithm 2.1) does not work when two or more final positions offer the same total cost for some given particle X i . In this case, the bidding increment γ i is zero, thereby creating a never ending cycle (see [5, Section 1.

3.3]).
To break such cycles, a perturbation mechanism is needed. Precisely, for a given small ε > 0, we will not look for fulfilling the conditions (2.3) any more, but the weaker ones, the so-called ε-complementary slackness, The Auction Algorithm is the same as the Naive Auction Algorithm, except that the bidding increment is 7) rather than γ i = −v i + w i as in (2.4).

Remark 2.2.
To make the connection with the assignment problem and the auction algorithm as they are presented in [6], we recall that we pass from our formulation to the classical one by taking a ij = −d ij .
• The cost formula (2.2) corresponds to the net value of object j for person i. Therefore, the increment γ i in (2.7) is equal to that in [6, Equation (8)].
The next result establishes the finite convergence of the Auction Algorithm (Algorithm 2.1 with (2.7) instead of (2.4)) to an approximate solution of problem (1.7). Proposition 2.3. The Auction Algorithm terminates in a finite number of iterations with a feasible assignment and a price vector that satisfy the ε-complementary slackness (2.6). Moreover, the overall cost of this assignment is within N ε of being optimal, i.e., if {(X 1 , Q j(1) ), . . . , (X N , Q j(N ) )} is the final assignment of the algorithm and D * = min (2.8) Proof. By Remark 2.2, Algorithm 2.1 with (2.7) instead of (2.4) takes the form of the auction process described in [6]. it follows from (2.8) that the assignment obtained at the end of the Auction Algorithm is optimal, i.e., it is the solution of problem (1.7).

An Improved Algorithm : the ε-Scaling
As explained in [5], one can improve the complexity of the Auction Algorithm by introducing the ε-scaling: rather than applying the algorithm for a fixed ε, one starts with a large value ε 0 and then successively reduce it up to some final ε K+1 , sufficiently small (like that in (2.9) for example). In addition, for each k ∈ {1, . . . , K}, successive auctions with ε = ε k are done for N − N k particles among the N such that (N k ) k is decreasing until 0.

Algorithm 2.5.
Step 0. All the particles are unassigned, the price vector is zero, ε 0 and θ > 1 are fixed. Let k = 1.
Step 1. Evaluate the number of particles to be assigned : Step 2. Use the Auction Algorithm with ε k = ε 0 θ k , until N − N k particles are assigned.
Step 3. Keep the price vector {p k j } 1≤j≤N −N k associated to the set of N − N k assignments among the N . These prices are used as input of iteration k + 1. While N k > 0 go to Step 1.
Setting θ > 1, let K > 0 be the greatest integer such that θ K ≤ N . Hence there exists q ∈ N such that q < θ K+1 − θ K and N = θ K + q. (2.10) Then, it is easy to see that Algorithm 2.5 stops at iteration K + 1 and that (the overall cost of) the final assignment is optimal within N ε K+1 . Therefore, the optimal solution is found for ε K+1 < 1 N in the case of integer costs d ij (see Remark 2.4). The case of real distances will be discussed in Section 3.
The main idea of the ε-scaling algorithm is to find an approximated solution for a coarse ε 0 (for which the prices increase rapidly at the beginning of the algorithm), and then to successively refine the solution with new (smaller) values of ε k . The final assignment, obtained with ε K , is exactly the same with or without considering the ε-scaling algorithm. But the claim is that the new algorithm should lower the total number of computational operations.  Integer K corresponds to the exponent in (2.10). Points in blue (K = 1) correspond to θ ≤ N < θ 2 , points in green (K = 2) correspond to θ 2 ≤ N < θ 3 .

Numerical experiments
Thanks to the forthcoming numerical experiments, we aim to study the influence of the following parameters, that have been introduced in Section 2: • ε which leads to the optimality condition (2.8). One can easily check that ε < 1/N is sufficient to obtain an optimal assignment in the case of integer distances d ij , but we need to investigate the case of real distances. • θ which is heuristic: is there a domain in R 2 of couples (θ, ε 0 ) for which the transport cost is minimum, leading to the solution of the assignment problem? Several criteria will be presented here through some experiments on integer and real data. Other properties that have to be considered are the complexity of the algorithm and its CPU cost (e.g. its computational cost on a workstation). For such an algorithm, which is strongly dependent on data configuration, we talk about the worse-case complexity: an upper bound for the complexity, obtained for bad-conditioned data. There is no theoretical study available so far, but it was evidenced (see e.g. [4,7]), that the worse-case complexity of the Auction Algorithm with ε-scaling and integer data is O (N 3 log(N C)), where Hereafter, we present the complexity study of our specific configuration, and then conclude on the efficiency and possible improvements of this algorithm.

The case of integer distances
The configurations of initial and final positions, called Test case 1, and shown in Figure 2, are the following: the initial particles are located at points (i, j, k) ∈ N 3 with 1 ≤ i, j, k ≤ 10, and the final one are translated by the vector (10, 0, 0). The assignment with minimal transport cost is then {((i, j, k), (i + 10, j, k))} i,j,k=1,2,..., 10 . In the configuration Test Case 1, we look for the couples (ε 0 , θ) leading to the optimal assignment. The results are given in Figure 3. This is done by running the Auction Algorithm on several number of particles N , and several couples (ε 0 , θ). If the numerical solution is exactly the optimal one, then a point is added in Figure 3. These points are classified with respect to colors, representing several values of K. The same N may appear several times, for instance N = 64 2 can also be written N = 8 4 . Several parameters may then possibly lead to the optimal solution, and the aim is to find an interval in which we are sure the optimal solution is always found.
As discussed in [5], a sufficient condition to reach the optimal assignment in the case of integer distances d ij = |X i − Q j | 2 is ε 0 < 1 (see Remark 2.4 and Equation (2.10)). This condition is not necessary, as can be seen in Figure 3.
According to [5], the choice of θ is heuristic (for integer distances as well as for real distances): 4 ≤ θ ≤ 10. For large N (e.g. for large K), the points are located in the interval 0 < θ ≤ 10, which means that for biggest θ, the algorithm did not find the optimal solution. In our tests, we thus found that the interval 4 ≤ θ ≤ 10 is not really well adapted: 4 ≤ θ ≤ 6 seems to be a more confident interval.

The case of real distances
The ε-scaling version of the Auction Algorithm aims only at accelerating the convergence of the algorithm, and does not change the optimality criterion of the solution. That is why we only deal in this part with the initial auction process for fixed ε (see Section 2.2 above).
This behavior is experimented in the configuration Test case 2, presented in Figure 4. The goal of this part is to find a link between N and the largest ε < 1 leading to the smallest transport cost. The interval 4 ≤ θ ≤ 6 previously deduced from the integer distances is also validated here, and in the following we take θ = 4.
Since we want to understand precisely the behavior of the algorithm in function of the disorder of the configuration, we only change the initial positions from Test case 1, taking initial random positions.
We first define an approximated optimal solution, using the following algorithm: (1) Choose a value for N .
(2) Consider a finite set of values of ε, compute the ε-optimal solution for each of them; we obtain a set of assignments S ε indexed by ε.
(3) For each assignment S ε , compute the associated transport cost D ε =    Figure 4). Each value of ε of Figure 5 is a mean of ten realizations for a given N . For a given N , most of the ε are equal, and the confident interval is O( ε). In this case, ε seems to follow a law like ε = O( 1 N ), which is the same behavior found in the integer case (see Remark 2.4).
Compared to figure 3, several remarks have to be made: first, in the case of integer distances, the optimal configuration is found even for high ε 0 , because of the simplicity of our test case. Second, the configuration of Test Case 2, which contains more disorder, implies lower values for ε. Finally, it is not sure at all that for all ε ′ 0 < ε 0 , the algorithm reaches the optimal solution. In figure 5, ε is the largest one found, ensuring that for all ε taken lower than it, the optimal configuration is found.

Complexity and CPU time
The worse-case complexity for the Auction Algorithm (with ε-scaling, in the case of integer data) is said to be in O(N 3 log(N C)) [4,7]. Surprisingly, Figure 6 shows a better result for Test case 2: the complexity in our case (according to ε, for a fixed θ), tends to be ruled by a law in O(N 2 ).
It is noteworthy that, in our application, the configuration and the regularity of the final positions of Test case 2 (shown in Figure 4) may improve the prefactor appearing in front of O(N 2 ). Indeed, taking a random repartition of the final positions rather than the regular one in Test Case 2 leads to the same behavior of the complexity, except that the prefactor is higher. The implementation of the Auction Algorithm uses efficient structures, but the time of execution in a sequential processor remains high: 3 min for 10 4 particles. The main limitation appears in the Step 2. of Algorithm 2.5: for large N , and a given number N − N k of particles to assign, the number of unassigned particles is not strictly decreasing, but can remain constant during numerous iterations of the loop 2.1. The reasons of such a phenomena are not well understood now, and we expect that some modifications (introduction of virtual source and well, for instance) may improve the Auction Algorithm.

Conclusion and forecast
After introducing the framework of the Stochastic Downscaling Method, we have presented an algorithm that corrects the particle positions. This algorithm, introduced by D.P. Bertsekas, is based on an auction process, and on an ε-scaling refinement; it contains several parameters. Our aim was to understand this algorithm, and use it as a reference method for the optimal transport problem.
We have proposed a way to characterize the optimality of the solution in the space of parameters, and a criterion on parameters ε 0 and θ which ensures to find a quasi-optimal solution in the case of real distances: 4 ≤ θ ≤ 6, and ε = O( 1 N ). It is noteworthy that, even if the algorithm does not necessarily provide the optimal solution, Equation (2.8) provides a control of the optimality error with the parameter ε.
There is no theoretical result available on the complexity of the algorithm, since its efficiency is strongly dependent on the data. However, in our case, the complexity is surprisingly low: O(N 2 ). Although the complexity behaves rather goodly, some work has to be done to bring down the CPU time, using improvements of the Auction Algorithm proposed in [4].
As a conclusion, this work validates the Auction Algorithm as a powerful tool for optimal transport problem in dimension three. This algorithm will be the reference method, that will be compared to other cheaper strategies solving the uniform density constraint, as the one proposed in [19]. Our future work thus consists in coupling SDM with the Auction Algorithm, and in adapting the Auction Algorithm to our application. Indeed, the resolution of the stochastic Lagrangian model is based on a Particle in Cell method, and this specific structure of cells should be integrated in the Auction Algorithm. This could lead to a cheaper solution in term of number of operations, and a more adapted solution, according to our physical requirements.