Thème Num

Marcelo Buffoni

Haysam Telib

Angelo Iollo

apport de recherche

Keywords: low-order models, domain decomposition, compressible flows modèles réduits, décomposition de domaine, écoulements compressibles

We propose a method to reduce the computational effort to solve a partial differential equation on a given domain. The main idea is to split the domain of interest in two subdomains, and to use different approximation methods in each of the two subdomains. In particular, in one subdomain we discretize the governing equations by a canonical scheme, whereas in the other one we solve a reduced order model of the original problem. Different approaches to couple the loworder model to the usual discretization are presented. The effectiveness of these approaches is tested on numerical examples pertinent to non-linear model problems including the Laplace equation with non-linear boundary conditions and the compressible Euler equations.

Introduction

In this contribution we are concerned with the coupling between a full order simulation and a reduced order model. The idea is to reduce the extent of the domain where we perform a canonical numerical simulation by introducing a low-order model which describes the solution far from the region of interest. By reducing the extent of the domain we aim at reducing the costs in terms of required memory as well as in terms of computational time.

In a broad sense, there exist many applications where far from the boundary the solution is weakly dependent on the details of the boundary geometry. In such regions we use a reduced order model based on proper orthogonal decomposition (POD) [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF] to solve the problem. This approach allows a representation of the solution by a small number of unknowns that are the coefficients of an appropriate Galerkin expansion. Therefore away from a narrow region close to the boundary of interest the number of unknowns to be solved for is drastically reduced. This idea was previously explored in the context of transonic flows with shocks [START_REF] Lucia | Reduced order modeling of a two-dimensional flow with moving shocks[END_REF], [START_REF] Legresley | Dynamic domain decomposition and error correction for reduced order models[END_REF]. Here we extend those works by adapting to that context some classical domain decomposition techniques.

In the following we discuss three possible methods to do that. The first is based on a Schur iteration where the solution of the low-order model is obtained by a projection step in the space spanned by the POD modes. The second is in the same spirit but instead of a Dirichlet-Neumann iteration we employ a Dirchlet-Dirchlet iteration in the frame of a classical Schwarz method. The last approach is of different nature since the solution of the low-order model is not simply based on a projection in the space of the POD modes. It takes into account in a weak sense the governing equations by minimizing the residual norm of the canonical approximation in the space spanned by the POD modes.

The numerical demonstrations shown in the following are relative to two models: the Laplace equation with non-linear boundary conditions modeling radiative heat transfer and the compressible Euler equations in a nozzle. Since this method can be of interest for optimal design applications, where many different geometries must be tested to improve performance, in some cases we have explored the idea of simulating by usual discretization methods just the region where the geometry changes, modeling the rest by POD.

Like all other approaches based on POD, a solution database is necessary to build the basis functions, therefore this method will be useful when many computations for relatively similar cases are to be performed, like for example in shape optimization, see [START_REF] Bui-Thanh | Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition[END_REF].

Solution by projection of the solution trace in the space spanned by the POD modes 2.1 Approximation of the Steklov-Poincaré operator by POD

In order to explain the method we take a particular case. Let us consider the Laplace equation ∆u = 0 defined inside a square

Ω = [0, 1] × [0, 1]. Let d ∈ [0, 1], Ω 1 = [0, d] × [0, 1], Ω 2 = [d, 1] × [0, 1]
and Γ = Ω 1 ∩Ω 2 the interface between the two sub-domains, see fig. 1. We have Dirichlet conditions on the right boundary (u R ) as well as on the upper (u U ) and lower (u D ) boundaries. Ω to Ω 2 . The coefficients b ik are found by POD as explained in [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. Parts I,II and III[END_REF]. It can be shown that this base gives by construction an optimal representation of the original data set u (k) Ω . Let us define û2 = ū2 + M i=1 a i φ i , where M is much smaller than the number of discretization points in Ω 2 . For an arbitrary Dirichlet condition on the left boundary of Ω, we want to determine the discrete solution by a canonical approximation in Ω 1 and by the above defined Galerkin representation in Ω 2 .

One simple way to implement this idea is to solve the problem by Dirichlet-Neumann iterations. To this end, we follow the steps below:

1. solve the problem in Ω 1 by any discretization method (FD, FEM, etc.), imposing Neumann b.c. on Γ;

2. on interface Γ project the trace of the above solution in the subspace spanned by the traces of the POD modes φ i ;

3. recover û2 as the prolongation of the trace of û2 on Γ inside Ω 2 by using the POD modes; This is just one possible solution algorithm, corresponding to a classical domain decomposition method (Schur complement). Another approach consists in solving the problem all at once, as detailed in the following. Let us define A 1 the discretized operator acting on u 1 , the restriction of the unknowns belonging to Ω 1 ; A Γ the discretized operator acting on u Γ , the unknowns belonging to Γ and A 2 the discretized operator acting on u 2 , the restriction of the unknowns belonging to Ω 2 .

The discretized non-linear problem in Ω can be written

  A 1 B 1 0 B t 1 A Γ B t 2 0 B 2 A 2     u 1 u Γ u 2   =   f 1 f Γ f 2   (1)
where B 1 and B 2 are appropriate interface matrices and f 1 , f Γ , f 2 take into account the boundary conditions.

From (1) we have

A 1 u 1 + B 1 u Γ = f 1 B t 1 u 1 + (A Γ -B t 2 A -1 2 B 2 )u Γ = f Γ -B t 2 A -1 2 f 2 (2)
The matrix A Γ -B t 2 A -1 2 B 2 is the discrete counterpart of the Steklov-Poincaré operator for Ω 2 , see [START_REF] Quarteroni | Domain Decomposition Methods for Partial Dierential Equations[END_REF]. Consider now the second step of the solution algorithm proposed above. Let a ∈ R M be a vector of components a 1 . . . a M and c ∈ R M a vector of components c 1 . . . c M . Posing ϕ i the trace of φ i on Γ, we take

a = arg min c∈R M u Γ - M k=1 c k ϕ k (3)
where ||•|| is the norm induced by the canonical l 2 scalar product, noted by (•, •). Solution of (3) reduces to the solution of the linear problem

M i=1 a i (ϕ i , ϕ j ) = (u Γ , ϕ j ), 1 ≤ j ≤ M . Therefore a i = (u Γ , P i ), where P i = M j=1 [(ϕ i , ϕ j )]
-1 ϕ j is a constant vector computed once for all from the POD modes.

At this point we approximate u 2 with û2 and substitute in [START_REF] Bui-Thanh | Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition[END_REF]. Since

B t 2 û2 = B t 2 ū2 + M i=1 a i B t 2 φ i , we have B t 2 û2 = B t 2 ū2 + M i=1 B t 2 φ i (u Γ , P i ). Finally, letting Ŝ2 = M i=1 B t 2 φ i P i we obtain the approximation of (2) B t 1 u 1 + (A Γ -Ŝ2 )u Γ = f Γ -B t 2 ū2 (4)
where

B t 2 ū2 ≡ B t 2 A -1 2 f 2 .
Matrix Ŝ2 is the approximation of the discrete Steklov-Poincaré operator obtained by the POD expansion. Equations ( 2) can of course be solved simultaneously by a standard linear solver.

Just like for the usual Steklov-Poincaré operator, (4) amounts to a non-local boundary condition for the problem posed in Ω 1 . The main advantage of this approach compared to computing explicitly S 2 is that we do not need A -1 2 to build Ŝ2 .

In the following we present some numerical applications. The first case is a plain application of the Schur complement as it was explained above to a non-linear case. The second case is based on the Schwarz method. The third numerical experience is a variant of the all at one method.

Schur complement

A second order finite differences (FD) method coupled to a fix point iteration is used to solve the Laplace equation inside the square domain shown in fig. 1. The left Dirichlet boundary condition is varied to build the needed database. In particular, u L = sin (kπy) + y, 1 ≤ k ≤ 49. The boundary conditions on the other sides are, referred to fig. 1, the following: on u U : u = 1, on u D : u = 0 and on u R : u 4 -u 4 0 + ∂u ∂n = 0. The domain is split at d = 1/3. Then, the POD basis functions are generated on Ω 2 using the previously computed database. In order to check the accuracy of the method, a boundary condition which was not included in the database used to build the POD modes is imposed on the left boundary: u L = y 2 , and a second order FD method is used to solve the problem in Ω 1 . We use 6 POD modes to recover û2 inside Ω 2 . Figure 2 presents the result of the test by means of the distribution of the relative error between the FD solution on the entire domain and the approximate Schur complement approach.

Schwarz method

In the following, a convergent-divergent domain Ω is considered. As before, Ω is divided in two subdomains, Ω 1 and Ω 2 , in a way that there exist an overlap region Ω ov = Ω 1 ∩ Ω 2 shared by both subdomains (see fig. For a geometry of the divergent part which is not included in the database, the solution is determined following a similar approach to that described for the Schur complement but employing this time the classical Schwarz method (see, for example, [START_REF] Quarteroni | Domain Decomposition Methods for Partial Dierential Equations[END_REF]):

1. solve the problem in Ω 1 by any discretization method imposing Dirichlet b.c. on Γ 1 n ; 2. on Γ 2 n project the trace of the above solution in the subspace spanned by the trace of the POD modes φ i ; 3. recover û2 as the prolongation of the trace of û2 on Γ 2 n inside Ω 2 by using the POD modes;

4. set u 1 = û2 on Γ 1 n ; 5. goto [START_REF] Bui-Thanh | Aerodynamic data reconstruction and inverse design using proper orthogonal decomposition[END_REF], n = n + 1, until convergence is attained.

Four POD modes are used to recover û2 inside Ω 2 . Figure 4 shows the results obtained for this case, again in terms of the relative error between a numerical solution (FEM P1) on the entire domain and the approximate Schwarz method.

Newton method

In this case the two-dimensional compressible Euler equations are solved

0 = -v • grad S 0 = -v • grad a -γ-1 2 a div v 0 = -grad v • v -a 2 γ-1 grad a -1 γ(γ-1) a grad S , (5)
with γ = 1.4 for air and being S the entropy, a the speed of sound and v = [u, v] the velocity, in a convergent-divergent nozzle on a structured mesh. We consider shockless flows. The λ-scheme [START_REF] Zannetti | About the numerical modelling of multidimensional unsteady compressible flow[END_REF] is used to solve the equations. Total temperature, total pressure and the flow angle are imposed at the inlet; static pressure at the exit and impermeability at the walls. The complete system is solved by Newton iterations. The resulting linear problems are solved by preconditioned GMRES iterations [START_REF] Saad | Sparsekit: a basic tool kit for sparse matrix computations[END_REF].

Using this code, a database of 90 snapshots is computed. The geometry of the divergent part of the nozzle is changed, Ω 1 in fig. 3. Furthermore, the static pressure at the exist is also varied taking uniformly spaced values in the interval p = [0.94, 0.99] with step 0.01. This corresponds to 15 snapshots for each pressure step. Since the flow is shockless, total temperature as well as total pressure are constant across the nozzle. Therefore, in order to solve in Ω 1 , the only unknown INRIA boundary condition is the flow angle on Γ 1 n . Hence, a low-order representation of the flow angle v/u for the convergent part, i.e., Ω 2 , is constructed retaining 10 POD modes.

We consider a case for which the geometry of the divergent part of the nozzle and the static pressure imposed at the exit do not belong to the set used to build POD database. The two subproblems are coupled by a least squares approximation on the overlapping domain. This is just a variant of the Schwarz method presented in the previous section. The only difference is that the two coupled problems are solved all at once by a Newton method.

The relative errors restricted to Ω 1 are reported in table 1 for some of the flow variables. These errors are computed with respect to the numerical simulation on the entire domain, and are quantified in terms of relative error in L 2 norm, i.e., the L 2 norm of the difference between the solution obtained by the present method and the reference solution divided by the L 2 norm of the reference solution.

Variable e(u)% e(v)% e(M a)% e(v/u)% Error 1.1 2.25 1.14 0.75 Table 1: Relative percentage errors (in L 2 norm) for the flow variables: u (horizontal velocity), v (vertical velocity), M a (Mach number) and v/u (flow angle).

Figure 5(a) presents the distribution of the Mach number obtained by the full numerical simulation on the entire domain. The same quantity obtained when applying the proposed method is shown in fig. 5(b). Finally, fig. 5(c) shows the distribution of the relative percentage error. In figure 6 the flow angle v/u is considered.

Solution by minimization of the residual norm in the space spanned by the POD modes

An alternative way to couple the low-order model to a detailed simulation is to look for an approximate solution in the reduced order function space that takes into account the governing equations. Hence, the main difference with respect to the approach described in section 2 is that the approximate solution is found by minimizing the residual norm of a given discretization scheme on the whole domain Ω 2 rather than by projecting the trace of the solution in the space spanned by POD modes.

Again the total pressure and the total temperature are constant across the nozzle and therefore we can write any other variable as a function of the local Mach number and the ratio v/u. In particular let U be the array of the couples (M a, v/u) for all the grid points belonging to Ω 2 . We start by representing this vector in the original Q-dimensional discrete space by a linear combination of basis functions, U = M i=1 α i Φ i , with M ≪ Q. The arrays Φ i have the same structure of U and they have been obtained by POD. The idea is to satisfy the compressible Euler equations (5) in a least squares sense over Ω 2 . In other words, let E(α) be the discrete residual of the governing equations as a function of α = (α 1 , . . . , α M ) and let be the residual norm. The solution in the POD function space α * is found by setting

I(α) = 1 2 E T (α)E(α) (6) RR n°6383 (a) (b) (c)
α * = arg min α I(α) (7)
This is equivalent to a system of non-linear equations

∂I ∂α i = ∂E T (α * ) ∂α i E(α * ) = J T (α * )E(α * ) = 0 ∀i ∈ [1, . . . , M ] (8)
where J T is the system's Jacobian. Equation ( 8) is solved by a quasi-Newton method where we take

E(α * ) = E(α 0 ) + J(α 0 )∆α + O(∆α 2 ) (9)
and after substituting in ( 8) we obtain

J T J∆α i = -J T E (10)
Therefore, we are left with the solution of equation ( 10) at each step of the quasi-Newton algorithm employed to solve [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. Parts I,II and III[END_REF].

Using the same database generated in Section 2.1.3, two low-order basis restricted to the convergent part of the nozzle are constructed. One for the flow angle v/u and the other for the Mach number. Only 5 POD modes are retained for both the flow angle and the Mach number.

This method is tested for the same case as for the Newton method presented in Section 2.1.3. The only difference is that the solution in Ω 2 is found by minimization of the residuals norm in the space spanned by the POD modes. Otherwise, the full numerical simulation is employed on Ω 1 , and a classical Schwarz overlapping method is used to iterate the solution to convergence.

The reference Mach number on the entire domain, the Mach number obtained by the present method on Ω 1 as well as the distribution of the relative percentage error on Ω 1 are shown in figures 8(a), 8(b) and 8(c), respectively.

In fig. 7, the results for the flow angle v/u are shown while table 2 These results show a much better accuracy as compared to those of section 2.1.3 as one can conclude by comparing tables 1 and 2.

Finally, in order to asses sensitivity of the results with respect to the minimization problem, we show in table 3 the error in Ω 2 when we use for α the initial guess. In table 4 we present the error at the end of the Schwarz iteration. It is seen that the fact of approximating the solution of the Euler equations with appropriate boundary conditions, even in a crude way, in Ω 2 , significantly improves the solution with respect to the initial guess.

Discussion

The fact of reducing the extent of the computational domain does not guarantee that one can get a solution faster as compared to solving the problem on the entire domain. In particular, let us consider the method described in section 2.1.3. It would in principle be the less demanding in terms of computational time since it just amounts to a non-local boundary condition in the frame of a Newton method. However, this is not necessarily the case since the convergence of the implicit iteration is spoiled by such non-local boundary condition and the number of Newton steps to attain convergence goes from 7 on the whole domain, to 11 when we only solve on Ω 1 . The slower convergence rate may lead to comparable costs in terms of CPU time, as a function of the fraction of the domain that is actually resolved. On the other hand, the fact of reducing the number of grid points is reflected almost proportionally on the memory requirements. However, the Jacobian matrix will have a nonsparse block corresponding to the non-local boundary condition induced by the approximation of the Steklov-Poincaré operator.

Concerning the method described in section 3 the cost of each step of the proposed minimization algorithm can be split up as follows. i) The computation of J. The Jacobian is evaluated by onesided finite differences. Consequently its cost in terms of CPU time and memory requirements is proportional to the dimension of the low-order space M and the cost of computing the residual vector, i.e., cost(J) = M × cost(E). Since M is O(10), the Jacobian matrix can be formed after few residual evaluations that are cheap in terms of CPU time. Memory requirements could become prohibitive for very big problems only on serial architectures. ii) The computation of the symmetric matrix A = J T J. The dimensions of A are M × M , with M = O(10). The required floating point operations are (M × M )/2 × N and they are less than the operations needed for the computation of E. iii) The cost of finding the solution to the linear system of size M × M can be neglected. Summarizing, it can be deduced that the cost is equivalent to some residual evaluations. The number of iterations needed in order to find the minimum are generally less than five, yielding an overall cost negligible with respect to canonical CFD calculation on Ω 2 .

We remark two major limits of these approaches. The first is of course that the results depend to a large extent on the database used for the POD modes. If the configuration under consideration lies in a region of the parameter space far from that explored when building the database, then the approximation error can be large. It is true, however, that the low-order model should be used where the solution does not strongly depend on the boundary conditions or the geometry. Another limitation is that an efficient way to improve the approximation quality comparable for example to grid refinement is not available. In principle, we would like to increase the approximation accuracy by enriching the functional space in which the solution is sought, based on some objective criteria. Unfortunately a general framework for such improvement is not presently available and it is the object of present research.

In conclusion we presented some possible implementations of a method to reduce the extent of the computational domain in the numerical solution of partial differential equations. The idea of using models that take into account different physical phenomena in different subdomains is old. Here we revisited this approach using a low-order model in the framework of classical domain decomposition techniques. The results in terms of the approximation error are promising for all of the cases that we showed. A major challenge to be pursued is to find a viable a posteriori error estimation technique for iteratively adapting the POD function space.

Figure 1 :

 1 Figure 1: Problem set-up: Schur

INRIA 4 .

 4 set ∂ û1 /∂n = ∂ û2 /∂n on Γ; 5. goto (1) until convergence is attained.

Figure 2 :

 2 Figure 2: Distribution of the relative error between the solution obtained by the present method and the solution obtained by a second order FD method on the whole domain. Ω 1 (a), Ω 2 (b).

INRIA

  3). The Laplace equation is solved with the boundary conditions detailed below to generate a database of k solutions with 1 ≤ k ≤ 60, this time varying the geometry of Ω 1 . The solution is obtained by means of the finite element method as implemented in [5], using P1 elements on a triangular non-structured mesh and a fixed point iteration. The boundary conditions imposed are: on u L : u = 1 3 • y, on u U : u = 1, on u D : u = 0 and on u R : u 4 -u 4 0 + ∂u ∂n = 0.

Figure 3 :

 3 Figure 3: Problem set-up: Schwarz

Figure 4 :

 4 Figure 4: Percentage error distribution on (a) Ω 1 and (b) Ω 2 . The reference solution is obtained by FEM P1 elements on a triangular mesh.

Figure 5 :Figure 6 :

 56 Figure 5: Distribution of M ach number: (a) Numerical solution (λ-scheme). (b) Present method. (c) Relative percentage error distribution on Ω 1 .

Figure 7 :Figure 8 :

 78 Figure 7: Distribution of v/u. (a) Numerical solution (λ-scheme) on the entire domain. (b) Present method. (c) Relative percentage error distribution on Ω 1 .

Table 2 :

 2 presents the relative errors obtained in Ω 1 . Relative percentage errors (in L 2 norm) for the flow variables: u (horizontal velocity), v (vertical velocity), M a (Mach number) and v/u (flow angle).

	Variable e(u)% e(v)% e(M a)% e(v/u)%
	Error	0.05	0.26	0.05	0.13

Table 3 :

 3 Relative percentage errors (in L 2 norm) for the flow variables: u (horizontal velocity), v (vertical velocity), M a (Mach number) and v/u (flow angle).

	Variable e(u)% e(v)% e(M a)% e(v/u)%
	Error	2.01	3.94	2.08	2.90

Table 4 :

 4 Relative percentage errors (in L 2 norm) for the flow variables: u (horizontal velocity), v (vertical velocity), M a (Mach number) and v/u (flow angle).

RR n°6383