C Vollaire
email: vollair@trotek.ec-lyon.fr

L Nicolas
email: nicolas@trotek.ec-lyon.fr

A Nicolas

PARALLEL COMPUTING FOR THE FINITE ELEMENT METHOD

A finite element method is presented to compute time harmonic microwave fields in three dimensional configurations. Nodal-based finite elements have been coupled with an absorbing boundary condition to solve open boundary problems. This paper describes how the modeling of large devices has been made possible using parallel computation. New algorithms are then proposed to implement this formulation on a cluster of workstations (10 DEC ALPHA 300X) and on a CRAY C98. Analysis of the computation efficiency is performed using simple problems. The electromagnetic scattering of a plane wave by a perfect electric conducting airplane is finally given as example.

INTRODUCTION

Nodal-based finite element (FE) method has been previously developed for microwave problems [START_REF] Nicolas | Three Dimensional FE Analysis of High Power Microwave Devices[END_REF]. Open boundary domains are modeled by coupling with Absorbing boundary conditions (ABC). The time harmonic formulation is written in terms of vector fields. Because nodal-based finite elements are used, a penalty term is added to the formulation in order to avoid spurious reflections. This code has been first developed on scalar workstation. This implies that only simple problems can be modeled. For instance, a 3 wavelengths side cubic geometry (30 cm at 3 GHz), meshed with 10 nodes per wavelength, leads to 81000 complex unknowns. With first order hexahedral FE, 181 Mbytes of memory are necessary to store the matrix. Obviously, the modeling of large geometry such as airplane can not be performed on scalar computers. Only parallel computation actually enables to modelize such devices: it reduces the computation time and, aboveall, it arranges enough memory.

The objective of this paper is to describe how the existing code has been modified in order to implement it efficiently on parallel computers. A cluster of ten DEC ALPHA workstations linked by a FDDI ring has been first used. This distributed memory computer is a multi instruction -multi data streams (MIMD) type. Parallel virtual machine (PVM) software is used to pass messages between workstations. The formulation has also been implemented on a CRAY C98 which is a shared memory computer (MIMD type) with vector capabilities.

Note that the algorithms and the methodology presented in this paper are not specific to a high frequency electromagnetic formulation. They can actually be applied to any physical problem discretized with a FE method.

We are dealing with frequency domain open boundary electromagnetic field problems. According to Maxwell equations, the magnetic field H and the electric field E satisfy to the vector wave equations. Following developments are made only for the H field formulation. All the steps can be applied to the E field formulation as well.

The weak Galerkin formulation of the vector wave equation for H is given by (1):

( ) ( ) ∇ × ∇× - ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ - × ∇× ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = ∫ ∫ W W v W s 1 1 0 0 2 ε μ ε r r v r s k d d H H n H (1)
A penalty term is added to the formulation to avoid spurious reflections [START_REF] Nicolas | Three Dimensional FE Analysis of High Power Microwave Devices[END_REF]. It makes the field divergence free. Its Galerkin form is given by (2):

[ ] [ ] -∇ ∇⋅ + ∇⋅ ∫ ∫ W v W s ( ) H n H v s d d (2)
ABC is used to truncate the 3D finite element region and to minimize the spurious reflections due to the outer boundary of the FE domain. A 3D vector Engquist-Majda condition [START_REF] Engquist | Absorbing Boundary Conditions for the Numerical Simulation of Waves[END_REF] is used. Outer surfaces are then rectangles, which is less mesh-consuming for a large number of geometries:

n H H H H × ∇ × ≅ = - ∇ g jk j k ABC t t t ( ) 0 0 2 2 (3)
Finally, the formulation for scattering problems in term of total field is given by (4).

H and H i are respectively the total field and the incident field. Two types of surface are considered:

-external surface (Sext ) -Perfect Electric Conductor (PEC) surface (Spec).

[ ] [ ] Moreover, the number of processors is limited by the decomposition method. This prevents the use of massively parallel computer.

( H ) H ( ) ( H ) n H (H) (H ) n H ∇ × ∇ × + ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ -∇ ∇⋅ + ∇⋅ - = -×∇× ∫ ∫ ∫ ∫ ∫ + W W v W v W s W s W s 1 0 2 ε μ r r v v Sext Spec ABC Sext Sext ABC i i k d d d g d g d (4)

IMPLEMENTATION ON THE CLUSTER OF WORKSTATIONS

Hence, we have preferred to perform the assembling step of the global matrix by degree of freedom [START_REF] Zois | Parallel Processing Techniques for FE Analysis: Stiffness, Loads and Stresses Evaluation[END_REF], [START_REF] Vollaire | Finite Elements Coupled with Absorbing Boundary Conditions on parallel Distributed Memory Computer[END_REF], instead of by elementary contribution as in a classical sequential code. It allows to assemble at once the three lines (three coordinates) corresponding to a node. Each element including the considerate node is sought and contributions with the related nodes are computed. These three lines are then compressed by storing only non-zero values. Space is left in the compressed storage if the volume element includes a surface element. Same method is used to introduce ABC on the external boundary and BC on conductors.

This method does not require the creation of the structure of the matrix (called symbolic assembling). To solve a problem meshed with N nodes, if P processors of the cluster are available, each workstation assembles a part of the system including N/P nodes. The load balancing is nearly perfect because of the constant bandwidth of the global matrix. The program is a SPMD type: no message passing is required. Hence, the speedup of the assembling stage is optimal (Fig. 3).

BC on symmetry planes and conductors

Messages passing are necessary to introduce the BC on conductors (only when solving for the E field) and on the symmetry planes because a global modification of the matrix system is required. This is due to the method used to introduce the BC (Fig. 1).

In the case of the H formulation the field verifies (6a) or (6b) on a symmetry plane:

On a symmetry plane: n H ⋅ = 0 (6a)

On an antisymmetry plane: n H × = 0 (6b)

The speedup related to this step is low due to the large amount of massage passing (Fig. 3). However, its effect on the global performances of the code is negligible: the length of this step is less than one percent of the total computation time.

Symmetrizing the global matrix

The introduction of ABC and penalty function leads to a non-symmetric system matrix. This one is approximately symmetrized by adding it to its transposed matrix.

This has been shown previously to give correct results. This operation requires a large amount of messages passing (Fig. 2). The speedup related to this operation is also low (Fig. 3). But once again, the time needed to perform this step is small and does not have any significant effect on the total computation time.

Solving in parallel

Because the FE matrix is sparse, iterative methods are used to solve the matrix system. Conjugate Gradient (CG) with diagonal, incomplete Cholesky or block incomplete Cholesky preconditioning methods have then been implemented [START_REF] Freund | Conjugate Gradient Type Methods for Linear Systems with Complex Symmetric Coefficient Matrices[END_REF], [START_REF] Joly | Mise en Oeuvre de la Méthode des Elements Finis[END_REF].

In parallel FE codes which use domain decomposition techniques to assemble the FE matrix, the solving step is first performed on the sub-domains and then on the global FE matrix [START_REF] Vollaire | Finite Elements Coupled with Absorbing Boundary Conditions on parallel Distributed Memory Computer[END_REF]- [START_REF] Lee | A Partitioning Technique for Finite Element Solution of Electromagnetic Scattering from Electrically Large Dielectric Cylinders[END_REF]. With our assembling method (degrees of freedom), the solver works on the entire FE matrix. So, we have chosen to use a small-scale parallelism.

Diagonal Preconditioning (DP)

Each processor computes his part of the preconditioning matrix by inversion of the diagonal terms of the FE matrix. No message passing is required and the memory space to store the preconditioning matrix is small. However, the generated preconditioning is low.

The preconditioning is performed by a multiplication vector-vector. A multiplication matrix-vector is also required to compute the residual vector at each iteration. This multiplication is done in parallel (fig. 4). Each processor computes a partial residual vector by multiplicating his part of the FE matrix by the vector duplicated on each processor. The load balancing is nearly perfect because of the constant bandwidth of the matrix.

To operate the concatenation of these partial residual vectors, each processor broadcasts it to all the others. Only non-zero values are sent. Then, each processor adds these partial residual vectors to obtain the final residue: this is the SPMD mode. This operation needs (P 2 -P) messages passing per iteration. An other way to calculate the residual vector consists in sending each partial residual vector to one processor called the master. This one (Master processor number 1) computes the final residue and broadcasts it. This Master-Slave (MS) method minimizes the number of communications: only (P-1)x2 messages per iterations are required. On the other hand, it introduces an additional idle time. Furthermore, because only non zero terms are sent, the messages broadcasted by the master are larger than in the SPMD mode: their size is equal to the size of the matrix. For these reasons the SPMD mode is more efficient. In any case, for both methods, the cost of communications is very penalizing, because small-scale parallelism is not adapted to distributed memory computer fitted with this kind of network. Some others methods can be found in [START_REF] Magnin | A Parallel and Vectorial Implementation of Basic Linear Algebra Subroutines in Iterative Solving of Large Sparse Linear Systems of Equations[END_REF] but they require more memory space.

Figure 5 shows an average of the state of every processor, on the total execution time, for both methods of concatenation (solving of a 60000 degrees of freedom problem on 4 processors). Figure 6 shows the corresponding speedups.

Incomplete Cholesky Preconditioning (ICP)

This preconditioning is performed by the decomposition of the FE in two matrices: A= L . L t (L and L t are computed). The building of the incomplete Cholesky matrix is performed by column (5) [START_REF] Magnin | Parallel and Vectorial Solving of FE Problems on a Shared-Memory Multiprocessor[END_REF], [START_REF] Magnin | A Parallel and Vectorial Implementation of Basic Linear Algebra Subroutines in Iterative Solving of Large Sparse Linear Systems of Equations[END_REF] (fig. 7). This algorithm is implicitly parallel because the L(ij) terms can be computed independently once the diagonal L(jj) term has been computed. The knowledge of both lines j and i is required to compute L(ij). If both lines are not stored on the same processor, a message passing is necessary. Fig. 7 and fig. 8 illustrate this strategy on an example of a 4x4 matrix stored on 3 processors. The incomplete Cholesky matrix is assembled in 5 parallel steps. Between each step, a processors synchronization is required.

For j = 1 to number of lines Ljj = Ajj -(Ljk)

For i = j + 1 to number of lines Lij = 1 Lii Aji - Lik Ljk End End 2 k=1 j-1 k=1 i-1 ∑ ∑ ( ) (5)
While solving using CG, both matrix-vector multiplication and back-forward substitution steps are necessary [START_REF] Zois | Parallel Processing Techniques for FE Analysis: Stiffness, Loads and Stresses Evaluation[END_REF]. They have also to be parallelized if possible.

The system A x = b is classically substituted by (L.L t ) x = b. Then, L y = b and (L t ) x = y are computed to solve the system. This algorithm is implicitly sequential because there is a back dependency on y and x: for example, the first processor computes his part of y and broadcasts it to all the others (in SPMD mode). The second processor can then start to compute and so on. This step is very penalizing in term of parallel performances, because it is performed at each iteration of the solver (fig. 9). Figure 10 shows the corresponding speedup.

Forward substitution

For i = 1 to number of lines (n) yi (bi Lik yk ) Lii k =1 k = i -1 = -∑ / End (6)

Back substitution

For i = number of lines (n

) to 1 xi (yi Lki xk ) Lii k = i +1 k = n = -∑ /

End

This method allows to reduce the number of iterations necessary to solve the system of equations (tab. 1). However, because of the large amount of message passing required, it cannot be applied to large problems. The memory space needed to achieve the ICP is 1.5 time larger than when using the DP because L and L t have the same structure than A after its symmetrization. L t is built to achieve a quick access by columns to the terms. This means that the terms are stored by lines in L and by columns in L t .

Block Incomplete Cholesky Preconditioning (BICP)

To avoid messages passing during the building of the incomplete Cholesky matrix, this one is assembled only with the terms stored on the considerate processor. So a part of the matrix is not built (fig. 11), and the terms effectively assembled are approximated because of the back dependency. This method applied with one processor corresponds to the classical ICP.

Compared to the Incomplete Cholesky Preconditioning, this scheme leads to a degradation (tab. 1) of the preconditioning. On the other hand, the preconditioning matrix is now constituted by independent sub-blocks. Each processor can compute independently his part of the result vector. The concatenation of these partial results, necessary at the end of the forward and the back substitution, is performed by message passing in SPMD mode. As show is tab. 1, the number of iterations required to solve the system of equations depends on the number of processors available: the increase of the number of processors leads to an increase of the number of iterations. On the other hand, the parallel rate of the method is increased (fig. 9-12). Figure 12 shows that the processors are better used than when solving with a classical ICP. The number of operations required to compute both preconditioning matrix and result vector depends on the number of processors available because of the not-computed terms in the preconditioning matrix. So, the parallel performances cannot be estimated in term of speedup because the solver on 2 processors is three times faster than on 1 processor.

Figure 13 shows the CPU time per processor for a solving of a 60000 degrees of freedom matrix versus the number of processors used.

For this size of matrix, the method loses efficiency when more than 5 processors are used. The increase of the matrix size tends to reject this limit.

Note that the memory space needed to achieve the BICP is 1.5 time larger than when using the DP because L and L t have the same structure than A after its symmetrization.

Comparison between methods

Table 1 compares the performances of the preconditioning methods when solving a 60000 degrees of freedom matrix on 2, 4 and 8 processors.

For this small example, on 2 processors, the Diagonal Preconditioning is the most efficient. The Block Incomplete Cholesky Preconditioning becomes interesting when more than 4 processors are used. For large problems, because of the large amount of messages passing, the building of the incomplete preconditioning matrix requires too much CPU time. Only the Block Incomplete Cholesky Preconditioning can be used.

Modeling a large problem

The electromagnetic scattering of a plane wave by a perfect electric conducting airplane (fig. 14) is presented as realistic problem. The frequency of the incident plane wave is 0.3 Ghz and this problem is meshed with 51183 nodes (307098 degrees of freedom). To solve the system of equations of this problem on 8 processors, the conjugate gradient with Diagonal Preconditioning requires 7476 iterations and 80280 s per processor while the conjugate gradient with the Block Incomplete Cholesky

Preconditioning requires only 4902 iterations and 58950 s per processor.

IMPLEMENTATION ON THE CRAY C98

The CRAY C98* is a high performance parallel vector computer. It is a shared -The first one consists in putting a parallelization directive before a loop. This one will be split, and the number of iterations performed by each processor can be specified. The variables used in the loop must have a scope, which means that they must be either private or shared. Each processor has its own copy of the private variable, while a shared variable is shared by all processors. Memory write conflicts on a shared variable may result in data corruption and must be avoided.

memory
-The second way consists in creating a parallel region, which includes a loop.

All the variables inside the parallel region must have a scope. The loop is parallelized, but the private variables can be used outside of the loop. This is useful for initializations, or to create private buffers to avoid memory conflicts. Semaphores can also prevent these conflicts.

A loop will be broken in packets that are long enough to enable vectorization (128 Words). Data streams must be as longer as possible to allow a good vectorization.

The CRAY compiler can try to parallelize and vectorize the code automatically due to compilation directives, by checking the data dependencies in the loops. It is automatic, but it works only on simple loops. For example, a subroutine call inside a loop prevents it from being parallelized. This automatic method leads then to bad parallel performances.

In order to parallelize loops which have not been made automatically, an higher level tool (BROUSE) may be used. It tries to find the scope of each variable in a specified loop. There are also mainly profiling tools as 'hpm', 'proview', 'atexpert', ... which give to the user many information about the performances of the code such as: parallelism ratio, average vector length, parallel / sequential portion of the code, most important loops / subroutines, ...

FE algorithm on the CRAY

The algorithms operating on mono-processor calculator do not need to be entirely revised in order to work on the CRAY C98. The code remains the same and runs quite immediately with more or less automatic parallelism-vectorization and progressively gains parallelism-vectorization while the user manually adds parallelizationvectorization directives.

Our code was first developed on a workstation, in a well structured way with many subroutine calls, error handling, .... Those points prevent the compiler from getting a good parallelization level. The compiler itself got only 17% of parallelization on a problem meshed with 10000 nodes. Only simple loops, like initializations had been processed. Neither the matrix assembly portion nor the solver had gained parallelism.

They only gained some vectorization. It was then necessary to help the compiler to get better results by adding manually compiler directives and by finding the scope of every variables.

The experience in the implementation of codes on such computer has shown that the matrix representation has a great influence on vector performances. So, two types of representations have been used.

A. Sparse row-wise matrix representation (storage # 1): because of the FE matrix is sparse and symmetric, only the non zero terms of its lower part are stored [START_REF] Engquist | Absorbing Boundary Conditions for the Numerical Simulation of Waves[END_REF] after its symmetrization.

B. Redundant Sparse row-wise matrix representation (storage # 2)

: because the algorithms used in the solver need to access to the FE terms by column, the entire FE matrix stays in memory even after its symmetrization. The memory space used is twice larger, but the access to the non zero terms of a column are adjacent in memory [START_REF] Magnin | A Parallel and Vectorial Implementation of Basic Linear Algebra Subroutines in Iterative Solving of Large Sparse Linear Systems of Equations[END_REF].

Some other methods like Sky line matrix representation have been implemented.

Because all the terms contained in the bandwidth of the FE matrix are stored, the number of terms stored for a line of the FE matrix increases with the size of the problem. This is due to the renumbering algorithm. So, this method requires too much memory and it can not be used to compute realistic devices.

All the performances are analyzed as previously with a 10000 nodes problem because the use of the ATEXPERT software (analysis of parallel performances) involves a mono-processor computation.

Creation of the global matrix structure

This step prepares the storage of the non-zero terms of the matrix. Its parallelization is difficult because there are many data dependencies. This part of code is especially vectorized.

Assembling by elementary contributions

This step is easy to parallelize. The elementary matrix related to a FE needs to have a private scope, and semaphores are used when the global matrix, stored in shared memory, is modified. So a parallel region is created which includes the loop on the FE.

The global matrix is modified in a critical region, inside the parallel region, in which the code is executed in a sequential way.

An other solution can be found in [START_REF] Ratnajeevan | FE Electromagnetic Field Computation on the Sequent Symmetry 81 Parallel Computer[END_REF] but it needs a 'pre-processing' step to split the set of FE into subsets made of independent elements. This method is efficient when write access conflicts in the global matrix for contiguous elements decrease the performances.

Introduction of BC on conductors (E field solving) and on the symmetry planes

Due to the method used to introduce BC, a global modification of the matrix system is required. Each processor performs a part of the modifications which are made on the global matrix. Semaphores are used to avoid memory conflicts.

Symmetrization of the matrix

Non-symmetric global matrix is symmetrized by adding it to its transposed matrix:

all the terms of each line are added with those corresponding of the transposed matrix.

The external loop (on the lines) is parallelized while the internal one (on the columns) is vectorized. The variables used for this step can be shared because no memory conflict is possible.

Figure 15 shows the speedups for the creation of the FE matrix (steps 4.

Solver

As previously presented, the conjugate gradient is used to solve the system of equations.

Diagonal Preconditioning

The use of this method requires only a multiplication matrix-vector per iteration to compute the new residual vector. So this step is parallelized by splitting the loop on the lines.

For the storage # 1, each processor computes a partial residual vector in private memory (parallel region). The addition of these partial results is performed in shared memory in a critical region (vectorization).

For the storage # 2, the multiplication is done in shared memory because the entire FE matrix is stored.

The speedups versus the matrix representation of the CG with the Diagonal Preconditioning are presented in figure 16.

Incomplete Cholesky Preconditioning

The storage # 1 requires the research of the lines (i) with a non zero term on the column (j). The parallelization is made by splitting this loop. Indeed, this algorithm needs the multiplication of the line (j) by the line (i). So an other search on the terms of the line (i) is necessary for all the terms of the line (j). This matrix representation is not adapted to the Incomplete Cholesky Preconditioning because of the vector performances obtained (tab. 2 and 3). On the other hand, only the lower part of the preconditioning matrix is built. Because only the lower part of the FE is stored in memory after its symmetrization, no more allocation of memory is required: the preconditioning matrix is stored in memory at the place of the upper FE matrix.

The storage # 2 allows to access directly to the lines (i) which have a non zero term on the column (j). Therefore, this loop is parallelized. The multiplication of the line (j) by the line (i) requires a search (for all the terms of the line (j)) on the terms of the line (i). This matrix representation improves the vector performances of the code during the building of the preconditioning matrix (tab. 2 and 3). The entire preconditioning matrix is built to keep the advantage of the adjacent access to the terms of a same column.

Because the entire FE matrix is stored in memory after its symmetrization, this type of representation requires to double the memory space.

Figure 17 shows the speedups obtained for the building of the Cholesky matrix versus the matrix representation.

The parallelization of the forward substitution cannot be made by splitting the loop on the line because of the back dependency. So, the produce of the terms L(ik) . y(k) is parallelized by creating a parallel region. The addition of the partial results is made in a critical region. The access to the terms L(ik) is performed by line, leading to good vector performances.

The parallelization of the back substitution is made in the same way but the access to the terms L(ki) decreases the vector performances (tab. 2 and 3). So the influence of the type of representation is the same as in the building of the Cholesky matrix. Figure 18 shows the speedups for back-forward substitution versus the matrix representation.

For the storage # 1, the research of the term L(ki) is penalizing in term of vector performances (tab. 2 and 3). The overhead introduced by splitting the loop on the produce L(ki) . x(k) is then negligible.

The use of the storage # 2 involves less efficiency in term of parallel performances because of the overhead introduced by splitting the produce L(ki) . x(k). On the other hand, because the terms L(ki) are adjacent in memory, the vectors performances are increased (tab. 2 and 3).

Comparison between the methods

Table 2 compares parallel and vector performances of the preconditioning methods when solving a 60000 degrees of freedom matrix with both types of matrix representation (total CPU time). Table 3 compares same performances of the preconditioning methods when solving a 307098 degrees of freedom matrix (scattering by a pec airplane). Both tables show the efficiency of the storage #2.

We have presented in this paper the implementation of an electromagnetic scattering code on parallel shared and distributed memory computers.

Cluster of stations

The parallel performances obtained when building the FE matrix are satisfying.

Because of the low preconditioning, the cost of communications with the Diagonal Preconditioning is penalizing. This is due to the small-scale parallelism which is not adapted to distributed memory computer equipped with this type of network. The Incomplete Cholesky Preconditioning allows to reduce the number of iterations but cannot be used on large devices because of the large amount of messages passing required. So, our new preconditioning method named 'Block Incomplete Cholesky Preconditioning' seems a good compromise in term of convergence rate and CPU time.

However, the memory space needed to achieve the BICP is 1.5 times larger than when using the DP.

Computed examples show that the modeling of realistic problem is still not possible on such a cluster because of the memory available. For example, the illumination of a fighter by a 3 GHz plane wave -high range radar freqency-would lead to a 10 6 nodes problem.

On an other hand, this code is immediately implementable on a computer such as the CRAY T3E. This parallel distributed memory computer is a MIMD type. Given its characteristics (256 processors and 32 Go of memory) it should allow to arrange enough memory to modelize large devices. Moreover, this machine is equipped with a very high performances network.

CRAY C98

The adaptation to a shared memory computer is easy because it does not ask a complete restructuring of program. The code remains the same and runs quite immediately. The CRAY compiler can try to parallelize and vectorize the code automatically due to compilation directives, by checking the data dependencies in the loops. It is automatic, but it leads to bad parallel performances. So the programmer has to manually add parallelization-vectorization directives and to cheek the scope of every variable to keep the control on the parallelism granulity.

The relative performances due to the vectorization are very low with a classical matrix representation method because the data streams are short (the matrix is sparse).

The use of the Redundant Sparse row-wise storage allows to obtain acceptable vector performances but it requires to double the memory space. Our experience shows that it is more important to favor the vectorization than the parallelization because every computation node of this machine is especially vectorial.

  multi instruction, multi data streams (MIMD) computer. It supports therefore a parallelism of strong scale. Indeed, every processor can execute a different program or a part of same program. All the processors have access to the global shared memory through a central connection. Main characteristics of the CRAY C98 are the following: Parallelism and vectorization are introduced by compiler directives. A vectorization directive leads the compiler to use the vector registers. There are several ways to use parallelism. Two of them are mainly used: * This work was supported in part by the Institut du Développement et des Ressources en Informatique Scientifique (CNRS).

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Message passing for the symmetry planes -Example of 4 processors.

  a32 x3 + a42 x4 a31 x1 + a32 x2 + a33 x3 + a43 x4 a41 x1 + a42 x2 + a43 x3 + a44 x4 x1 x2 x3 x4 a11 x1 + a21 x2 + a31 x3 + a41 x4 a21 x1 + a22 x2 + a32 x3 + a42 x4 a31 x1 + a32 x2 + a33 x3 + a43 x4 a41 x1 + a42 x2 + a43 x3 +

Figure 4 .Figure 5 .Figure 6 .Figure 7 .

 4567 Figure 4. Parallel matrix-vector multiplication-Example of a 4x4 matrix on 2 processors.

Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .Figure 13 .

 8910111213 Figure 8. Building of L and L t in 5 parallel steps.

Figure 14 .

 14 Figure 14. Perfect electric conductor aeroplane enlighted by a plane wave (magnitude of the magnetic field H).

Figure 15 .Figure 16 .Figure 17 .Figure 18 .

 15161718 Figure 15. Speedups for the assembling step.

Table 1 .

 1 CPU time per processor and number of iterations versus the preconditioning method and the number of processors.

	Number of processors	2		4		8
		time (s) ite. time (s) ite. time (s) ite
	Diagonal Precond.	415	204	367	204	316	204
	Cholesky Precond.	1781	59	901	59	1773	59
	Block. Cho. Precond.	1133	59	353	61	234	81
			CPU time (s)	Mflops	Iterations
	Diag. Precond. storage # 1	263		20		232
	storage # 2	187		40.2	232
	Cho. Precond. storage # 1	5314		2.3		76
	storage # 2	733		10		76

Table 2 .

 2 Total CPU time, vector performances and number of iterations versus the preconditioning method and the matrix representation (solving on 8 processors).

	CPU time (s) Mflops	Iterations

Table 3 .

 3 Total CPU time, vector performances and number of iterations versus the preconditioning method and the matrix representation (solving on 8 processors).


Olivier Bourgeois
email: olivier.bourgeois@neel.cnrs.fr

Dimitri Tainoff

Adib Tavakoli

Yanqing Liu

Christophe Blanc

Mustapha Boukhari

André Barski

Emmanuel Hadji

Reduction of phonon mean free path: from low temperature physics to room temperature applications in thermoelectricity

Keywords: phonon, thermal transport, thermoelectricity, semiconductor, nanowire, inclusions phonon, transport thermique, thermoélectricité, semiconducteur, nanofil, inclusions

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.

Introduction

In the two last decade, materials structured at the nanoscale having low thermal conductivity have attracted a lot of attention due to their high potential thermoelectric applications [1,2,3,4]. Indeed, as proposed by Dresselhaus et al. in 1993, nanostructuration can significantly improve the thermoelectric efficiency [5,6]. This efficiency of thermoelectric materials is given by the dimensionless figure of merit ZT = S 2 T σ/k where S is the Seebeck coefficient, σ the electrical conductivity and k the thermal conductivity. If we focus only on the thermal properties, the ZT will be increased by reducing the phonon thermal transport. This reflects the fact that the heat source and the cold source must be as thermally isolated as possible to prevent too much heat link that will degrade the overall performances of the thermoelectric modules.

Since concerning electron, the thermal and the electronic properties are coupled through the Wiedemann-Franz law, the only way of reducing the thermal transport requires the reduction of the phonon thermal conductivity. Two major directions are available to act on the phonon properties: incoherent or coherent effects. On one hand, incoherent effects are ruled by the phonon mean free path (MFP) then the reduction is obtained by adding scattering centers, by reducing the size, having rough surface or by adding nano-inclusions. On the other hand coherent effects will manifest themselves through the reduction of the group velocity (band flattening), opening of band gaps or having destructive interferences. Two significant length scales govern the phonon physics involving incoherent or coherent effect: the phonon mean free path (Λ) and the dominant phonon wavelength given by:

λ dom = hvs 2.82kB T (1)
where h the Planck constant, kB the Boltzmann constant and vs is the speed of sound. This length is indeed temperature dependent, increasing as the temperature decreases (see Fig. 1).

Here, we will focus our attention on the experimental attempts to reduce of phonon mean free path in single crystal semiconducting materials by playing on the geometry of the surface or on the internal structure. Other means based on coherent effects will be discussed at the end of the paper. There is actually two complementary approaches to realize an efficient phonon engineering at the nanoscale: the top-down one based on e-beam lithography to make small structures like nanowire and membrane and the second one, the bottom-up, generally based on specific growth of nanostuctured materials (bulk or nanowire) using for instance molecular beam epitaxy (MBE).

These two important limits will be experimentally addressed in this paper summarizing the recent findings of our group: the low temperature one set by boundary scattering limit where both important lengths (Λ and λ dom ) play a role and the room temperature limit where for the time being only the mean free path can be engineered efficiently.

Towards the ballistic limit: thermal transport in Si nanowires

In low dimensional systems, the surface scattering of phonons plays a crucial role in the understanding of thermal transport [7,[START_REF] Ziman | Electrons and Phonons[END_REF][START_REF] Berman | [END_REF]. The so-called Casimir regime, occurs when the bulk phonon MFP is much bigger than the actual size of the nanosystem under study. In that case, the phonon MFP is only set by the scattering on the surfaces. This limit is reached at low enough temperature when other scattering mechanisms do not dominate (phonon-phonon interaction, electron-phonon interaction, scattering on impurities etc...). As a result, effects linked to surface roughness have a significant impact on the phonon transport [10,11]. These effects can be simple incoherent scattering when the phonon wavelength is much bigger than the mean asperity or, as the temperature is lowered, may become wavelength dependent. Indeed, if the wavelength is bigger than the mean roughness, specular reflections of phonons may occur on the surface. If elastic scattering is involved, the energy is conserved leading to an increase of the MFP. Phonon spectral density as a function of the wavelength for different equilibrium temperature of a black body. The dominant phonon wavelength for each temperature corresponds to the maximum of the curves. Typically the dominant phonon wave length is of the order of 1 nm at 300 K and 100 nm around 1 K.

Normal phonon transport in Si nanowires at low temperature

In order to understand the thermal transport at low temperature a summary of the Casimir and Ziman models is necessary. In the Casimir model, [7,[START_REF] Berman | [END_REF], the thermal conductance is given by:

KCas = 3.2 × 10 3 2π 2 k 4 B 5 3 v 3 s (2/3) e × wΛCas L T 3 = βCasT 3 , (2)
where kB is the Boltzmann constant, the Planck constant and vs refers to the sound velocity. ΛCas = 1.12 √ e × w is the Casimir MFP of the phonons in a nanowire having a rectangular section, e refers to the thickness and w to the average width of the nanowire; L being its length. This formula is obtained when the mean free path of phonons is only limited by the cross-section of the wire, i.e. by boundary scattering. Indeed, in this limit, each phonon hitting the rough surface is assumed to be absorbed an re-emitted in all direction and energy following the black body radiation law (see Fig. 1). In other words, in this Casimir limit, the surface is considered as infinitely rough whatever the phonon wavelength.

As the temperature is lowered, the dominant phonon wavelength actually grows (see equation 1) and becomes eventually larger than the roughness amplitude; then specular reflections become more probable and cannot be neglected like in the Casimir model. To take this trend into account, the Casimir mean free path should be replaced in equation 2 by an effective MFP Λ ef f as proposed by Ziman and coll. [START_REF] Ziman | Electrons and Phonons[END_REF][START_REF] Berman | [END_REF]:

Λ ef f = 1 + p 1 -p ΛCas ( 3
)
where p is the parameter describing the probability for a phonon to be specularly reflected at the surfaces. The value of this parameter will depend on both the wavelength of the phonons and the roughness of the wire surface. The Casimir purely diffusive case is given by p = 0 and the purely specular one, when p = 1, will be referred to as the ballistic limit; the MFP is diverging to infinite. The Ziman model takes into account a value for p that is different from zero, as the temperature is lowered the phonon mean free path can be bigger than the section of the nanowire. The different phonon transport regimes and their associated models are summarized in the Table given in Fig. 2. As a conclusion, the thermal transport in silicon nanowires having section of 100 nm will be well described by the Casimir and then at lower temperature by the Ziman model as it will be seen in the next section. Indeed, the length L is smaller than the bulk phonon mean free path Λ ph-bulk in Si, and on the same order of magnitude as the dominant phonon wavelength. At even lower temperature, the fully ballistic limit should be described by a Landauer type model where the thermal transport is given by transmission coefficients between the nanowire and the heat bath. This goes clearly beyond our focus in this paper; there is currently no experimental data in this regime.

Phonon blocking with S-shaped chicanes in Si nanowires

We will illustrate in the following section these different scatterings and phonon transport regimes in three situations: firstly the thermal transport in a normal straight silicon nanowire, secondly in a Si nanowire having Sshaped chicanes, and thirdly in Si nanowires with engineered surfaces. The nanowire thermal conductance has been measured using longitudinal 3ω technique using a highly sensitive thermometry based on niobium nitride (NbN) [12,10,11,13]. The Si nanowires are fabricated from silicon on insulator (SOI) substrates by e-beam lithography. The total length of each structure has been purposely set to 10µm in order to easily compare the thermal measurement of all kinds of nanowires (see Fig. 3). The section of the nanowires is 100 nm by 200 nm. The central symmetry of the structures is necessary for the 3ω method employed to measure the thermal conductance [12].

The temperature profile looks like a parabola with the hottest point located at the center of the nanowire. The heat generated at each points of the suspended structure flows from the center to the heat bath on both sides. The thermal conductance is generally measured using a power of a few tens of femtowatts dissipated in the niobium nitride transducer deposited on top of the nanowires. This power creates on average a temperature gradient smaller than 1 mK. The values of the measured thermal conductance are around a few universal value of thermal conductance K0 = π 2 k 2 B T /3h ∼ 10 -12 T W/K 2 ; a value used to normalized the data. The equation 2 gives the correct order of magnitude for the thermal conductance since a K close to 10 -11 W/K is predicted at 3K and a value of 1.5 × 10 -11 W/K is measured on a straight nanowire. For the explanation of temperature variation of the thermal conductance it is needed to take into account the contribution of specular reflection by introducing the Ziman mean free path, and hence the parameter p. This has been described in details in the past [10,11,14]. We would like to highlight here that the contribution of ballistic phonons to the thermal transport is significant. We have proposed a solution for thermal management that will not perturb diffusive transport but will affect only ballistic one. Specific suspended nanowires have been made including a double bend chicane structures of 400 nm long; the total extended length of the S-shaped nanowire itself being identical to the one of the straight nanowire (see Fig. 3c and Fig. 4).

In recent works [14,15], we have shown that the S-shaped nanowires have a much lower thermal conductance than straight one, as it can be seen in Fig. 4. If the transport would be only diffusive, this constriction would have no effect at all on the thermal conductance. Actually, the severe reduction observed (reduction by more than 30 %) demonstrates two significant facts: first, this illustrates the contribution of ballistic phonon to the thermal transport and second, including a chicane in a thermal conductor is a very efficient way of reducing the thermal transport. Phonons having long mean free path (bigger than the Casimir mean free path ΛCas) will be blocked and reflected by the chicane acting like a mirror for the long wave length acoustic plan waves. Inversely, short wavelength phonons that are more scattered by the nanowire surfaces won't be much affected by the chicane; their mean free path being already set by the section of the nanowire like in the Casimir regime. This concept has been successfully used to thermally isolate a suspended platform dedicated to sensitive thermal measurements [16].

Reduced thermal transport in corrugated nanowires

In the previous section, we have seen how a macroscopic constriction can act as a phonon filter. In this section, we will illustrate how phonon backscattering or multiple scatterings can be introduced by using nano-engineered surfaces. In this respect, phononic crystals, defined as a regular modulation of the geometry of a phonon waveguide, have attracted significant attention in the past decade as a promising means to reduce phonon transport. In the present samples, periodic nanostructurations are used (see Fig. 5) to induce phonon backscattering and reduce incoherently the thermal conductance; not coherently like in a real thermocrystal. We have reported that in that kind of engineered corrugated nanowires the phonon transport can be decrease even below the Casimir limit [15,17].

In order to compare thermal transport between various geometries of nanowires, the thermal conductances have to be normalized to a geometrical factor called the Casimir coefficient βCas, according to Eqs. 2 and 3. Fig. 5 reports the thermal conductance plotted as K/βCas. This normalization yields the temperature power law and allows a direct comparison between the thermal transport of each nanowire. Indeed, from Eqs. 2 and 3, the normalized conductance can be expressed as: where

log(K/βCas) = a + b log(T ) (4)
(K/βCas) = 1 + p 1 -p T b (5)
From these two latter equations, we can extract the power law factor b of the temperature variation of the thermal conductance, the value of the MFP Λ ef f and the parameter p [15,17].

The thermal conductance of the corrugated nanowires is clearly below the thermal conductance of straight nanowires even if the power law has not been change; being very close to 3 (b ≈ 2.6). The presence of artificial periodic nanostructuration increases the scattering of phonon (backscattering, phonon trapping etc...). The surface roughness is not changed but the mean free path is clearly reduced. This can be depicted by ray-tracing Monte Carlo numerical simulations. Indeed, the multiple scattering of ballistic phonons can be seen through a corpuscular picture explaining the fact that the phonons are somehow trapped in the core of the nanowire. As a consequence, the phonon MFP that can be even smaller than the Casimir limit, i.e. smaller than the mean diameter of the nanowire. This is characterized by the necessity of putting a negative p parameter to fit the data. Even if this special result does not involve phonon band engineering like in a real phononic crystal, it gives however a fantastic way of turning a single crystal semiconductor into a phonon glass.

In the diffusive limit: nano-inclusions limits thermal conductivity in semiconductor

After having shown in details the different ways of affecting the thermal transport in the presence of ballistic phonons, we will focus now on room temperature effects where the diffusive regime is governing the thermal conductivity. We will show that engineering the phonon thermal properties is also possible with semiconductor at room temperature (RT). As compared to the ballistic limits, working at RT in the diffusive limit will have significant implications for real applications in managing heating, cooling or harvesting energy.

The great challenges mostly rely on our capacity of improving the semiconductor materials performances. As on one hand a semiconductor is by far too thermally conductive (around 100 W/mK), on the other hand its Seebeck Figure 6. TEM image of the germanium matrix with GeMn spherical nano-inclusions at 8 % of Mn. The presence of numerous inclusions will significantly increase the phonon scattering without perturbing the electronic transport. The scale bar represents 50 nm. The distance between the nano-inlcusions is indeed smaller than the average phonon mean free path.

coefficient is generally quite high (above 200 µV/K). Then obtaining an efficient thermoelectric material having a high ZT based on a semiconductor is possible by decreasing significantly its thermal conductivity.

The increase of the thermoelectric efficiency is feasible thanks to the difference of mean free path existing between phonons (∼100 nm at 300 K) and electrons (∼1 nm at 300 K). By introducing in a given material a structural disorder at the nanometer scale, it is possible to induce phonon diffusion without affecting the charge carriers justifying the innovative concept of electron crystal-phonon glass material.

Indeed, in a semiconducting matrix, phonons contribute the most to the heat transport having an average mean free path of ∼100 nm at room temperature and a dominant wavelength of 1 nm at RT. With the objective of manipulating and controlling the phonon transport in a SC, the engineering of materials can then follow the two different approaches mentioned in the introduction: either by playing with phonon wavelength [18,19,20,21] or with the phonon mean free path [17,20,21]. At room temperature, the most efficient effect will be the reduction of the phonon mean free path using a geometrical scattering approach. Indeed, as it has been recently experimentally demonstrated the phonon wavelength is far two small at RT for the nanostructuration to create sizable effects [20,21].

The best large scale applicability is expected from bulk 3D SC material that has been engineered at the low dimension for thermal management. One of the most important concepts is the electron crystal-phonon glass being an inclusion-embedded semiconducting system, which shows great advantages and potentials [23]. For crystalline semiconductors, the ideal model would be a defect-free highly crystalline matrix containing nanoscaled inclusions with a wide diameter distribution. It can be theoretically predicted that as the inter-distance of the inclusions is comparable to the phonon mean free path in the matrix, the nano-inclusions can prohibit efficiently the phonon transport. Indeed, by introducing more phonon scattering processes we expect to get a germanium based electron crystal-phonon glass material.

The materials of interest here are thin films of a single-crystalline germanium matrix embedded with randomly distributed Ge3Mn5 nano-inclusions grown by epitaxy on a germanium substrate (see Fig. 6). The nano-inclusions are found to be nearly spherical with a diameter distribution of from 5 to 50 nm, depending on the growth parameters. Furthermore, manganese atoms act as p-type dopant in the Ge matrix, which ensures a good electrical conductivity of the material.

The experimental measurements of the thermal conductivities of the Ge:Mn materials have been done using a highly sensitive differential 3ω technique [22]. Experimental results of 3ω measurements revealed strong reductions of the thermal conductivity, compared to the Ge bulk at room temperature, for Ge:Mn thin films containing different Mn concentrations.

The major experimental finding is a reduction of more than a factor of ten from the bulk value in the 8 % of Mn in the Ge matrix. A thermal conductivity below 8 W/mK has been measured which has to be compared to the bulk thermal conductivity 65 W/mK. The reduction is fully ascribed to the presence of evenly distributed phonon scatterers than significantly reduce the phonon mean free path.

Conclusions

Different concepts of phonon thermal transport reduction have been evidenced experimentally from the ballistic limit to the diffusive limit. These new paradigms will be fruitfully applied for the development of new materials or devices in thermoelectricity and energy harvesting in general. It seems that now, we are facing a kind of intrinsic limit in that reduction by non-coherent effects, a lower limit close to the alloy value of thermal conductivity. Unfortunately, engineering the phonon band by opening gaps like in a phononic crystal seems not perfectly appropriate for room temperature thermal management, phonon wavelength being far too small as compared to the length accessible by nowadays fabrication means.

Other technological solutions may be found in playing also on the electrical properties of semiconductor based systems to improve the Seebeck coefficient. Indeed as proposed by Hicks and Dresselhaus [5] it could be possible to increase the Seebeck coefficient by playing on the electronic density of state at the Fermi level. Then, engineering also the electronic properties by various means (nanostructuration, adding quantum dots, etc...) seems absolutely needed if one wants to boost the performances of future thermoelectric materials.

Figure 1 .

 1 Figure1. Phonon spectral density as a function of the wavelength for different equilibrium temperature of a black body. The dominant phonon wavelength for each temperature corresponds to the maximum of the curves. Typically the dominant phonon wave length is of the order of 1 nm at 300 K and 100 nm around 1 K.

Figure 2 .

 2 Figure 2. Classification of the different regimes of phonon transport in nanoscaled systems at low enough temperature, typically below 30K for silicon. L is the sample length, Λ ph is the inelastic phonon mean free path, and λ dom the dominant phonon wave length. Concerning the Landauer limit, the parameter T denotes the transmission coefficient between a nanowire and its heat reservoir. Other important lengths are explained in the text.

Figure 3 .

 3 Figure 3. a SEM images of: a) a suspended nanowire linked to the two pads where the electrical connections permit the 3ω measurements b) a straight nanowire c) corresponds to the side view of a S-shaped nanowire.

Figure 4 .

 4 Figure 4. Thermal conductance of 5µm long nanowires normalized to four times the universal value of thermal conductance versus temperature. The thermal conductance of the straight nanowire is much bigger than the thermal conductance of the S-shape nanowires. The presence of the serpentine section stops the transmission of ballistic phonons, reducing the overall phonon mean free path.

Figure 5 .

 5 Figure5. Normalized thermal conductance versus temperature for a straight nanowire (Straight) and four corrugated nanowires (C1 to C4) in a log-log plot. The normalization coefficient βCas corresponds to the Casimir limit value of the thermal conductance at 1 K. The corrugated nanowire Ci can have much lower thermal conductance that smooth nanowire. The phonon MFP is reduced severely, even below the Casimir limit. As an illustration, in insert we compare the average phonon mean free path in the straight nanowire (in yellow), in the Casimir limit (in red), and as measured in corrugtaed nanowire (in blue). The corrugations reduce significantly the phonon mean free path.

Acknowledgments

We acknowledge technical supports from Nanofab, the Cryogenic and the Electronic facilities and the Pole Capteur Thermométrique et Calorimétrie of Institut Néel for these experiments. Funding for this project was provided by a grant from La Région Rhône-Alpes (Cible), by the Agence Nationale de la Recherche (ANR) through the project QNM no. 0404 01, by the European projects: MicroKelvin FP7 low temperature infrastructure Grant no. 228464 and MERGING Grant no. 309150.


Alain Jean-Marie

Mabel Tidball

Extraction cost: before or after harvesting? Economic and environmental consequences

Introduction

The model Solution Qualitative properties

Context of the paper

Groundwater extraction: marginal cost depends on the level of the aquifer In general, resources with accessibility problems: cost depends on scarcity

The main ingredient: make the cost depend on the projected evolution of the resource: before or after the extraction or rainfall

The goal: deduce its economical and environmental consequences The method: revisit the discrete time, innite horizon dynamic game model of Provencher and Burt (JEEM 1993).

Introduction

The model Solution Qualitative properties

Contributions

We analyze a variant of Provencher and Burt's model, in the linear-quadratic case.

We characterize the existence of Nash equilibria in ane feedback We prove qualitative properties of equilibria as a function of the discount factor.

We discuss the particular cases: myopic (zero discount factor) and golden rule (discount factor tending to one).

We focus on the case of scarse resources and nd that taking harvesting and rainfall into account in the cost is better than the standard situation.

The model of Provencher and Burt

We consider the extraction of groundwater by two players.

The dynamic of groundwater:

G t+1 = G t + R u 1 t u 2 t ; G 0 ; given:
We suppose R is a constant.

The instantaneous prot:

i = F i (u i t ) C i (G t ) ¢ u i t :
The marginal extraction cost (C i (:)) depends on the current level of the groundwater. The Bellman equation associated with a Nash Feedback equilibrium is:

V i (G 0 ) = max u i t ¢ F i (u i ) C i (G + mR n(u 1 t + u 2 t ))u i + V i (G t + R (u 1 t + u 2 t ))
£ :

F i (u) = u b i 2 u 2 ; C i (x ) = z i c i x > 0:
We propose as solution of the Bellman equation: a value function for Player i P f1; 2g of the form

V i (G ) = A i 2 G 2 + B i G + C i ;
a feedback law of the form u i = i G + i :

The unknown A i , B i , C i , i , i , i = 1; 2, are found identifying the coecients of the quadratic function in the Bellman equation after optimization the optimal strategy for j and the one conjectured by i. 9 / 29 The model Solution Qualitative properties General linear-quadratic solution The symmetric case Solution in the Linear-Quadratic case (ctd) More precisely, a solution by stages: system of 3rd degree polynomial equations for f 1 ; 2 g A 1 and A 2 as simple functions of f 1 ; 2 g

A i = c i (1 j n) i (b + 2c i n) (1 1 2 )
linear system for fB 1 ; B 2 ; 1 ; 2 g simple formulas for fC 1 ; C 2 g.

Interior solutions

The existence of a useful solution is not granted because:

the LQ problem is not concave i (u) = u(1 z i b i 2 u) + c i uG
G t = (1 1 2 )G t 1 + R 1 2 with solution: G t = (1 1 2 ) t G 0 + R 1 2 1 + 2 (1 (1 1 2 ) t ):
Necessary conditions for the trajectory to be valid are:

0 1 + 2 < 2 1 + 2 R.
Otherwise, an interior solution of our problem does not exit.

Introduction

The model Solution Qualitative properties General linear-quadratic solution The symmetric case

The symmetric linear quadratic case Symmetric case: look for symmetric equilibria: 1 = 2 = .

With the procedure described above we nd that solves: p(Z) := ¨2(b + 2cn)Z 3 (3b + 8cn)Z 2 + 2cZ © +( 1) f (b + 3cn)Z + cg = 0:

Existence and uniqueness of solution

There is one unique root of p(Z) in (0; The green golden rule, = 1 Note that the limit when goes to 1 can allow to select a solution for the static green golden rule max u i F i (u i ) C i (:)u i ; such that u 1 + u 2 = R:

A = c (b + 3cn) (1 2) = c 1 z + mcR R[c (b + 3cn)]
As it is well know, this is a game with coupled constraints and there exists an innite number of solutions.

In the symmetric case, lim t3I u i t = R=2 for i = 1; 2 and all .

But in the asymmetric case this limit can allow to select one equilibrium.

Observe that lim 41 (1) T = 0. Benchmark situation: the myopic case = 0 When = 0, we nd:

u(G) = c b + 3cn | {z } (0) G + 1 z + cmR b + 3cn | {z } (0)
:

The value function of both players is:

(G 0 ) = (cG 0 + 1 z + cmR) 2 (b + 3cn) 2 b + 2cn 2 :
And the asymptotic stock is:

G I = R 2(0) 2(0) = R(b + 3cm 2cn) 2(1 z) 2c : Introduction The model Solution Qualitative properties
Variation with respect to the discount factor Comparison of cost situations The myopic case The dynamic case Clearly: @(0) @n < 0; @(G 0 ) @n < 0; @(0) @m > 0; @(G 0 ) @m > 0:

Accordingly:

Ranking of controls, gains and steady states Controls are ordered as: u H0 (G ) < u 00 (G ) and u H0 (G ) < u HR (G ); u 00 (G ) < u 0R (G )

Value functions are ordered as: H0 (G ) < 00 (G ) and H0 (G ) < HR (G ); 00 (G ) < 0R (G )

Stochastic extension

Assume the recharge is a i.i.d. sequence fR t ; t = 0; 1; : : :g.

In the LQ case, the Bellman equation becomes:

V i (G 0 ) = max u i t F i (u i ) C i (G

  not depend on z, R or m and B depend linearly on m and R

  u H0 (G ) < u HR (G ) < u 00 (G ) < u 0R (G ) H0 (G ) < HR (G ) < 00 (G ) < 0R (order to have positive harvesting policy for all G and a positive steady state, we must impose after harvesting and rainfall is a good economic and environmental option Assume: = 0:3; b = 1; c = 0:04; z = 0:9; R = 1: interest to charge users in function of their behavior, not just in function of the level of resource Possibility of win-win situations More analysis to better explain the

The model Solution Qualitative properties Variation with respect to the discount factor Comparison of cost situations The myopic case The dynamic case

Qualitative analysis

We investigate the variation of the equilibrium (feedback parameters, value) when the discount factor varies when the cost adjustment parameters n (harvest) and m (rainfall) vary.

Introduction

The model Solution Qualitative properties Variation with respect to the discount factor Comparison of cost situations The myopic case The dynamic case Variation with respect to When agents are more shortsighted, they react more agressively and the environment suers, whatever the cost structure.

Monotony in

Under the assumption c < b=2, the function:

We conjecture that (¡) is actually always decreasing.

Introduction

The model Solution Qualitative properties Variation with respect to the discount factor Comparison of cost situations The myopic case The dynamic case

Variation of and

Towards a generalization of these ndings to general P [0; 1].

Variations of and

The function n U 3 (n) is decreasing, so that, 0R = 00 > H0 = HR :

The function m U 3 (m) is increasing, so that: H0 < HR ; 00 < 0R We have: when > 0, and R small enough, the total ranking of controls hold: u H0 (G ) < u HR (G ) < u 00 (G ) < u 0R (G ):

when $ 0 and z < 1, then > 0: the equilibrium control u(G) is valid for all G when $ 1, < 0. More analysis is needed for large and R small. u 00 (G ) > u HR (G ); V 00 (G ) > V HR (G ); when G big, same ranking as in the myopic case: conict; u 00 (G ) < u HR (G ); V 00 (G ) < V HR (G ); when G small, reversed ranking: win-win.

When the level of the groundwater is small, setting costs after harvesting and rainfall is better from the economic and environmental point of view than the standard literature case where the cost is announced before rain and harvesting, even for quite myopic agents.