Yves Mah

Jean-Louis Pazat
email: pazat@irisa.fr

Distributed Array Management for HPF Compilers

High Performance Fortran and other similar languages have been designed as a means to express portable data parallel programs for distributed memory machines. As data distribution is a key-feature for exploiting the parallelism of applications, a crucial point for HPF compilers is their ability to manage e ciently distributed arrays. We p r e s e n t in this paper an innovative method to allocate local blocks and temporaries for received values and to manage the associated access mechanisms. The performance of these access mechanisms is measured and experimental results on the use of this array management within an existing compiler are shown.

Introduction

In order to alleviate the task of programming Distributed Memory Parallel Computers, new features have been added to sequential programming languages such a s F ortran. In the eld of scienti c programming, two main axes are currently followed. The rst one uses explicit parallel constructs and loop partitioning, it relies on a shared virtual memory 2] the second one is based on a user-de ned data distributions which are used as a guideline to generate communicating processes 7].

In recent y ears, many projects have focused on the data distribution approach and it has been demonstrated that \aggressive" optimizing compilers and e cient r u n time systems are mandatory to achieve reasonable speedups. Most compilers allow the user to specify a decomposition of arrays and use the owner write rule 4] to distribute the code. This distribution can be done using the runtime resolution technique in which e a c h statement o f the source program is guarded and where communication is performed elementwise. This scheme is always applicable but rather ine cient, so that many compilers integrate optimization techniques for compiling loops. Roughly speaking, these techniques aim at reducing iteration domains and performing vectorized communications 14, 1 2 , 1 1 , 1 0 ] .

However, runtime resolution as well as optimization techniques require a speci c and e cient runtime system to allocate local parts of arrays and to perform e cient communications.

In this paper we present a new method for e ciently managing distributed arrays (allocation and access) in parallelizing compilers based on data distribution. This paper focuses on the local management o f b l o c ks and temporary storage for distant v alues. Communication optimizations such as message coalescing and vectorization are supported by this array management but are not addressed here.

The paper is organized as follows: next section discusses the essential requirements for managing distributed arrays. Section 3 details the page-driven array management proposed. Section 4 presents in more details the implementation of this array management w h o s e overall performance is presented in section 5. Future work is discussed in the conclusion.

2 Management of Distributed Arrays

Key issues

The aim of a distributed array management is to de ne the way distributed arrays are stored in the local memories and the way elements of these arrays are accessed. These tasks involve both the compiler and the runtime support. It is clear that a trade-o must be achieved between the speed of accesses and the memory overhead induced from the array representation. The extreme solution consisting in allocating the entire array o n e a c h processor is obviously not applicable. Conversely, m i n i m al allocation (typically achieved by translating an array declaration A(N) i n to a local allocation A(N=P) where P is the number of processors) associated with global-to-local index conversion involving several costly operations such a s mod and div must be avoided.

In addition to the \classical" concerns that are the memory use and the speed of accesses, we de ne in the following some properties that we claim to be useful for a distributed array management.

Uniformity

Two kinds of data are to be considered on a given processor: local elements { i.e. elements assigned to it by the distribution { and received elements that are temporarily stored in the local memory after communication with another processor. We s a y that an array management s c heme is uniform when the representation as well as the access mechanism associated with a distributed array m a k es no distinction between local and received elements. The main advantage of a uniform scheme is that the compiler does not have to separate purely local computation (involving only local data) from computation that needs distant d a t a . Indeed, even if such a separation is sometimes possible at compile-time, it may induce an important code fragmentation in the case of multiple right-hand-side references performing the separation at runtime brings about costly ownership computation.

Furthermore, de ning a non-uniform scheme that gives too great an importance to local elements, at the expense of received elements is not a good solution because it is likely to be harmful not only for accesses but also for communication of non-local data.

Independence

An independent array management is only de ned from distribution parameters and in particular does not depend on the code itself nor on compiler analysis on the code. A distributed array management s c heme independent from the compilation scheme facilitates the coexistence of di erent compilation techniques. On the contrary, i f t h e c hoice of a representation is guided by the analysis of a program part (typically a loop), it may h a p p e n t h a t s e v eral layouts (and associated access methods) are used within the scope of one distribution, possibly necessitating data rearrangement or additional computation at runtime. Moreover, this independence facilitates the use of di erent compilation techniques within a code fragment that contains several loops. One way t o a c hieve this independence (as well as uniformity) is to consider that only global indices appear in the generated code.

Contiguity Preservation

Another useful property concerning the layout of distributed arrays is the conservation of memory contiguity. Indeed, if contiguous elements of the original array are still contiguous in the local representation, it makes it possible to take a d v antage of direct communications (in this case no copying nor packing/unpacking is needed between local representations and communication bu ers), vector processors, target code optimization and better cache behavior.

Related Work

To our knowledge, management of distributed arrays have not been studied independently from compilation techniques in existing HPF compilers.

The rst technique of storage for distributed arrays, the overlap 13] has been implemented both in the Vienna Fortran Compilation System 14] and in the Fortran D compiler 8, 12]: a single sub-array is allocated for local data as well as for received data. This technique provides uniform accesses and preserves memory contiguity but it can be applied to a restricted number of distributions and access patterns and may lead prohibitive allocation when distant data location is not close to the local partition. In this case, the Fortran D compiler may select an alternative storage method (bu ers) for received values if it can separate purely local computation and computation needing received values. In Vienna Fortran and in an extended version of Fortran D, a speci c management related to loops with irregular array accesses is performed through the use of the inspector/executor technique: the local partition is dynamically extended during the inspector phase so that uniform indirect accesses can be used during the executor phase [START_REF] Chatterjee | Generating Local Adresses and Communication Sets for Data-Parallel Program[END_REF][START_REF] Bodin | Fortran-S : A Fortran Interface for Shared Virtual Memory Architectures[END_REF][START_REF] Zima | SUPERB: A Tool for Semi-Automatic MIMD/SIMD Parallelization[END_REF].

Other array management s c hemes, closely related to compilation techniques, have been proposed. Their common characteristic is that they try to minimize memory overhead. Among them, in the compilation scheme de ned by Ancourt et al 9], local elements and temporaries are packed according to the array distribution and alignment, the loop bounds and the array subscripts by c hanging the basis of the original index space. Accesses to elements are performed in a non-uniform way with index conversion evaluating a ne functions and possibly integer division. Chatterjee et al. 5] propose an access mechanism for local elements based on a Finite State Machine (FSM). These elements are accessed by executing a FSM that has to be computed at runtime for each l o o p n e s t e v en if the same distribution applies.

All of these methods not only take i n to account the array distribution parameters but necessitate also a static analysis of code fragments (loop bounds and array subscripts) in order to de ne the layouts of the local arrays and the associated access mechanisms. Therefore, independence from the compilation scheme is not achieved by these systems.

A P age-driven Array Management

We present here a new management s c heme for distributed arrays based on software paging. This management is designed in order to achieve e cient accesses while avoiding unacceptable memory overhead. It also aims at satisfying the properties of uniformity, independence and contiguity preservation aforementioned. In the following, we will consider only direct HPF distributions, the mechanisms described can be easily extended to aligned distributions (as a rst step, by applying paging to templates).

Principle

The page-driven data management w e propose follows the main addressing scheme of classic paging systems for memory management. In such systems, logical memory space is broken into groups of contiguous elements (pages). Pages have a xed predetermined size. A hardware support divides a logical address in two parts: a page number and a page o set. The page number is used as an index into a page table that contains the base address of each page in physical memory. This base address is combined with the page o set to de ne the physical memory address. If the page size is S, a logical address produces a page number PG and an o set OF by PG = d i v S and OF = m o d S . If the logical address space is larger than the physical address space, virtual memory management features may be added. In this case, accessed pages may not be present all the time in physical memory but temporarily loaded from a secondary storage system by a s w apping device.

As for our concern, we manage variables |i.e. distributed arrays| and not memory our aim is not to build a shared virtual memory. Moreover, we only consider compiler-generated code, hence we stay at the software level rather than relying on hardware components. Contrarily to system-level paging, the notion of page fault is here irrelevant because all distant accesses are solved by prior communications. Besides, data are not necessarily communicated page-wise. The original address space is multidimensional therefore we apply a multi to onedimensional transformation before splitting the resulting space into pages. This allows us to de ne a speci c access mechanism for each distributed array, in particular the page size may be di erent for each a r r a y.

Paging Distributed Arrays

We de ne a representation and its access mechanism for each distributed array b y a c o u p l e (L S ). The multidimensional index space of a given array is linearized by a function L. T h e linear address space obtained is split into pages of xed size S. A processor stores only those pages that contain at least one element a s s i g n e d t o i t b y the distribution or one received element. Depending on the distribution of the array, L and S, a page may be possessed by one or several processors.

One of the main advantages of this method is that accesses to local and received elements are performed the same way. Indeed, as far as accesses are concerned, a processor acts as if the entire array w as allocated locally, no matter if the element it needs to access is truly local or has been received from another processor. The di erence between pages containing local elements and pages containing received elements lies in the way they are allocated and lled, not in the way they are accessed. A tuple (PG,OF ) is computed from the initial index vector (i 0 : : : i n;1 ) with the linearization function L and the page size S: PG = L(i 0 : : : i n;1 ) div S OF = L(i 0 : : : i n;1 ) mod S A table of pages is stored on each processor. It indicates the base address of each p a g e present in local memory. The o set is added to this base address to obtain the exact location of the element. The page partitioning is also used for computing owners of elements. A similar table, present in the local memory of each processor, stores for each page the numbers of the physical processors that own this page.

Tuning Parameters

For a given distributed array, the parameters we can tune for paging are the page size S and the linearization function L. The value of these parameters are de ned in order to achieve good performance in terms of time and memory space.

As speed of access is our prior motivation, time consuming operations (division, modulo and multiplication) are avoided in the computation of the tuple (PG,OF ) but also in the application of function L. T h i s i s a c hieved by i n troducing powers of two, turning integer division, modulo and multiplication into simple logical operations. Moreover, the array d ecomposition can be taken into account when xing the actual value of S and L. I n tuitively, we c hoose S and L so that pages \follow" the blocks, and are owned by as few processors as possible.

For a more formal de nition, let us consider the following HPF array distribution: REAL V(0 : h 0 ; 1 0 : h n;1 ; 1) !HPF$ DISTRIBUTE V(CYCLIC(s 0 ),: : : ,CYCLIC(s n;1 )) and the access to an element o f V noted V(i 0 i n;1 ). This distribution decomposes the array V into rectangular blocks of size s 0 s n;1 . W e consider the most general distribution directive CYCLIC(k). Note that BLOCK(k) and CYCLIC(k) distributions are strictly equivalent as far as decomposition is concerned.

Prior to the de nition of S and L, w e c hoose a particular dimension , the dimension in which the block size is the largest. If there are several such dimensions, the one corresponding to a non-distributed array dimension or a block size equal to a power of two i s c hosen. The page size S is then given by: if s = h or s = 2 then S = sup (s ) else S = inf (s ) where sup (x) (resp. inf (x)) extends an integer to the smallest (resp. largest) power of two greater (resp. less) than or equal to x.

L is the C linearization function for multidimensional arrays applied to a permutation of the index vector. This permutation puts the index corresponding to dimension in last position. Moreover, the array dimensions (coe cients of L) are extended to the next power of two. L is de ned by L(i 0 : :

: i n;1 ) = n;1 X k=0 i 0 k n;1 Y l=k+1 h 0 l !
where i 0 k is the k th access index after permutation, i.e: i 0 n;1 = i 8k 2 0 : : : ;1 i 0 k = i k 8k 2 : : : n ;2 i 0 k = i k+1 and h 0 k is the extended size of the array i n t h e k th dimension, i.e:

h 0 n;1 = sup h S S if n > 1 h 0 0 = h 0 if > 0 else h 1 8k 2 1 : : : ;1 h 0 k = sup (h k ) 8k 2 : : : n ;2 h 0 k = sup (h k+1 )
Here are two examples of de nition of S and L rst when there is one non-distributed dimension and second when all the dimensions are distributed: L(i j) = 512j + i

Optimizing the Computation of (PG,OF)

Unlike with a classic paging mechanism, the explicit computation of the linear address L(i 0 : : : i n;1 ) before its splitting into (PG,OF) is not mandatory because we do not rely on a hardware support that needs a memory address. Besides, this intermediate result may lead to unnecessary operations as in the following example:

REAL A(0:99, 0:199) !HPF$ DISTRIBUTE A(CYCLIC(10), ) ) S = 2 5 6 L(i j) = 256i + j

The page number and the o set will be obtained by PG = (256i + j) div 256 OF = (256i + j) mod 256 These expressions could obviously be simpli ed in PG = i and OF = j. To m a k e t h e simpli cations clearly visible, we express directly PG and OF as a function of the index vector.

page(i 0 : : : i n;1 ) = ( PG OF) with PG = n;2 X k=0 i 0 k n;1 Y l=k+1 np 0 l ! + i 0 n;1 div S OF = i 0 n;1 mod S where np 0 k is the number of pages in the k th dimension after permutation: np 0 n;1 = h 0 n;1 S 8k 2 0 : : : n ;2 np 0 k = h 0 k When dimension is not distributed, that is to say w h e n h = s , index i 0 n;1 (i.e i ) i s always less than or equal to S, div and mod can be removed: Each processor stores a table of owners that indicates, for each p a g e , t h e n umber of the physical processor that owns this page. This table can be lled using the function owner(PG,OF ) that returns the owner of an element. owner(PG OF) = map page ;1 (PG OF) Function page ;1 , the reverse function of page, returns the index vector corresponding to a page number and an o set.

page ;1 (PG OF) = ( i 0 : : : i n;1 ) with i = S ; PG mod np 0 n;1 + OF 8k 2 0 : : : ;1 i k = i 0 k 8k 2 +1 : : : n ;1 i k = i 0 k;1 8k 2 0 : : : n ;

2 i 0 k = PG mod n;1 Y l=k np 0 l ! div n;1 Y l=k+1 np 0 l !!
Function map(i 0 : : : i n;1 ) associates a physical processor number with an index vector. It can be easily computed for each mapping that can be expressed in HPF. We do not give here the general formulae as it cannot be induced from the HPF norm and depends on implementation choices. As an example if the abstract processor array i s o f s i z e ( p 0 p n;1 ) and the number of physical processors is P, the mapping function may be the following:

map(i 0 : : : i n;1 ) = n;1 X k=0 (i k div s k ) n;1 Y l=k+1 p l !! mod P
In the case data elements are replicated (by application of alignment directives), this function could return a set of processors. The de nitions adopted for S and L allow t h e n umb e r o f o wners of a page to be less than or equal to two. If the owner of a page is always unique, any v alid value of OF can be used for determining the owner of a page. In the case the owner of a page is not unique, we can compute OF lm , the o set from which the owner changes. In this case, the table of owners stores for each page, the two processor numbers plus the limit OF lm :

OF lm = i f ' < S then ' else 0 with ' = ( ( PG mod np 0 n;1 ) (s ; S)) mod s

Implementation

A full implementation of the data management m e c hanisms described above has been realized within the Pandore environment 1 ]. In the Pandore language, the scope of array distributions is con ned within procedures referred to as distributed phases.

Management of tables, pages and accesses to array e l e m e n ts are shared out among the compiler and the runtime library. As all the tables and pages are needed only during the execution of a distributed phase (no inter-phase analysis is performed at this time), the entire memory space allocated is freed at the end of the phase.

Tables and Pages

All the information needed to ll the tables of owners and the tables of o set-limits is known at compile-time these tables could therefore be statically de ned. However, in order not to lengthen the size of the generated code, the compiler produces functions that allocate and ll the tables at runtime, at the beginning of each distributed phase. For each distributed array V , a table of owners TO V is de ned. If a page may be possessed by t wo processors, three tables are needed: the table of the owners of the rst part of pages TO1 V, the table of the owners of the second part of pages TO2 V and the table containing the o set-limits TL V.

The runtime library is also in charge of allocating and lling the tables of pages and pages themselves. The tables of pages and pages that contain local elements are allocated at the beginning of the distributed phase. The management of pages containing received elements depends on the compilation scheme. Basic operations provided by the runtime library are the page allocation and the placement of elements (single elements or segments) into pages.

Accesses

It is clear that the part of the access process that is done at compile-time m ust be as large as possible. The compiler translates a reference to an array element V I], where I is an index vector, into a call to a runtime macro access(desc V, PG, OF) where PG and OF are expressions of I. All constant subexpressions have been computed and the optimization described in section 3.4 has been performed. As expected, these expressions contain only additions and constant logical shifts and maskings. The work that remains at runtime is therefore to evaluate the expressions and use the table of pages associated with V (TP V) t o produce the right reference. This can be noted by the C expression *(TP V PG]+OF). T h e runtime library contains cpp macros that prevent from the computation of the address of the page table corresponding to V , s o w e can actually generate this code.

Owner

Determining the owner of an element V I] is carried out a similar way. The compiler generates a call to a runtime macro owner(desc V, PG, OF). An access to a table TO V PG] is su cient a t r u n time to nd the processor number in the case the owner of a page is unique. If a page may be possessed by t wo processors, a call to a slightly di erent m a c r o is produced. The execution of this macro will issue a comparison between OF and the o set-limit corresponding to page PG: if (OF < TL_V PG]) then TO1_V PG] else TO2_V PG] 5 Performances

Performances of the Distributed Array Management

It is quite obvious that the executed code for distributed accesses involves only few basic operations that generate a very small overhead and may e v en be more e cient thanks to better optimizations. In the experiment whose results are reported in table 1, we measured the time taken by several kinds of read accesses: t s : a reference to an element a s i t m a y appear in a sequential program t p : a call to the macro that uses the paged access mechanism t b : a call to a macro that uses a block-oriented access mechanism1 . The array i s a t wo-dimensional array of oats reported times are in s. Best and worst cases have been considered, depending on whether the sizes of the array w ere powers of two or not. Experiments have been carried out on a SparcStation 2, on a node of the iPSC/2 and on a node of the Intel Paragon XP/S. Likewise, the determination of the owner of an array element requires only a few simple operations, so its cost remains very low. It is also preferable to exploit the page decomposition, although it seems to be more natural to base the computation of the owner of an element on the computation of the corresponding block n umber.

The price to pay for speed of access and speed of ownership computation is the need for a larger amount of memory. O v erhead is only due to tables because no additional space is required for pages. When a page contains elements that will never be accessed because of the extension of array dimensions, or because the page is shared by t wo processors, only the potentially accessed part of the page is actually allocated. A translation of the corresponding pointer in the table of pages is performed if the end of the page is allocated.

The memory overhead due to tables is directly linked to the number of pages, which is in general at least of an order of magnitude less than the size of the array. Table 2 gives memory requirements for a few common distributions of arrays on 32 processors. For each distribution, we indicate the total number of pages, the theoretical minimal memory space required on each processor, the actual space allocated for tables on each processor and nally the overhead as compared with the minimal partition. Memory needs are expressed in bytes. It can be noticed that replacing some block sizes (or array dimensions) by p o wers of two notably decreases the memory overhead. We believe t h a t o verall memory requirements remain acceptable when considering most commonly used distributions.

Integration in the Pandore Environment

The page-driven management for distributed arrays has been integrated in the Pandore environment and is used with the two compilation schemes of the compiler. The basic compilation scheme, that relies on a runtime resolution technique, can be applied to every input program. The optimized scheme is based on integer programming and linear algebra results it performs an analysis of parallel loops 10].

We p r e s e n t in table 3 the results of the execution of a Red-Black Successive O v er-Relaxation algorithm run on a 1024x1024 matrix of oats. Times have been measured on the Intel iPSC/2 for the two compilation schemes. A comparison is made between a blockoriented array management and the page-driven management. The table shows the speedup obtained on P processors for each pair (compilation scheme, array management). The use of the page-driven management clearly improves performances of codes generated according to both compilation schemes. The joint use of the optimized scheme and the page-driven array management leads to satisfactory performances (e ciency of around 75% for 32 processors) in spite of the unfavorably high ratio of memory operations to computation of the Red-Black algorithm.

Conclusion

Management of distributed arrays is a crucial point for obtaining good performances when using data parallel compilers. For this purpose, we h a ve proposed a new scheme based on parameterized software paging that proved e cient in an existing compiler. This management handles local and received data in an uniform way and it is independent f r o m t h e optimization techniques used in compilers. Moreover, it avoids using multiple representations of the same array in di erent parts of a program and maintains some regularity i n local layouts. The page-driven array management also seems to be appropriate for irregular computations and could be used together with the inspector/executor technique.

We are currently comparing our management s c heme with shared virtual memory systems. One of our objectives is to nd out which features should be added to existing shared virtual memory systems so they can be e ciently targeted by data-parallel compilers.

  is the result of these optimizations for the two examples presented in the previous section:REAL A(0:199, 0:99, 0:49) !HPF$ DISTRIBUTE A(CYCLIC(5), , CYCLIC(10)) ) 8192i + 128k + j) div 128 = 64i + k OF = (8192i + 128k + j) mod 128

Table 1 :

 1 Speed of page-driven access

	Sparcstation best worst best worst best worst iPSC/2 Paragon ts 0.30 0.42 0.94 2.05 0.16 0.26 tp 0.34 0.38 2.14 2.26 0.22 0.25 tb 0.48 1.58 3.52 9.86 0.21 2.68
	Array Distribution REAL(KIND=8) A(0:99999) !HPF$ DISTRIBUTE A(CYCLIC(1000)) REAL(KIND=8) A(0:99999) !HPF$ DISTRIBUTE A(CYCLIC(1024)) REAL(KIND=8) A(0:999, 0:999) !HPF$ DISTRIBUTE A(CYCLIC(1), ) REAL(KIND=8) A(0:999, 0:1999) !HPF$ DISTRIBUTE A(CYCLIC(50), CYCLIC(500)) REAL(KIND=8) A(0:999, 0:1999) !HPF$ DISTRIBUTE A(CYCLIC(50), CYCLIC(512)) REAL(KIND=8) A(0:99, 0:99, 0:99) !HPF$ DISTRIBUTE A( , CYCLIC(1), CYCLIC(50))	Number Minimal Local Space of Pages Partition for Tables Overhead Local 196 25000 1960 1:08 98 25000 588 1:02 1000 250000 6000 1:02 8000 500000 80000 1:16 4000 500000 24000 1:05 10000 250000 60000 1:24

Table 2 :

 2 Memory overhead for some common distributions

Table 3 :

 3 Comparison between block-oriented and page-driven managements

	Basic scheme Optimized Scheme P Block Page Block Page
	4 0.29 0.68 0.86 8 0.32 0.77 1.39 16 0.36 0.82 2.14 32 0.37 0.85 3.70	3.84 7.18 12.78 23.72

This mechanism was used in a previous version of Pandore it performs at runtime a modulo and an integer division to nd the block n umber and the o set in the block.


Franw Andr6

Marc Le Fur

Yves Mah6o

Jean-Louis Pazat

The Pandore Data-Parallel Compiler and its Portable Runtime

This paper presents an environment for programming distributed memory computers using High Performance Fortran. Emphasis is put on compilation techniques and distributed array management. Results are shown for some well known numerical algorithms. 1

Introduction

The difficulty of programming massively parallel architectures with distributed memory is a severe impediment to the use of these parallel machines. In the past few years, we have witnessed a substantial effort on the part of researchers to define parallel programming paradigms adapted to Distributed Memory Parallel Computers (DMPCs).

Among then J, the Data Parallel model is an interesting approach: the programmer is provided a familiar uniform logical address space and a sequentiM flow of control. He controls the distributed aspect of the computation by specifying the data distribution over the local memories of the processors. The compiler generates code according to the SPMD model and the links between the code execution and the data distribution is enforced by the owner-writes rule: each processor executes only the statements that. modify the data assigned to it by the distribution. This approach constitutes the basis of several compilers [START_REF] Tseng | An Optimizing Fortran D Compiler for MIMD Distributed-Memory Machines[END_REF][START_REF] Zima | Compiling for Distributed-Memory Systems[END_REF] and is also applied in the PANDORE compiler.

This paper presents the PANDORE environment and focuses on some optimizations of the compilation scheme and the run-time support. It is organized as follows: the next section presents the PANDORE environment. Section 3 explains the compilation process and the optimizations used for parallel nested loops. The array management is described in section 4 and results are discussed in section 5. Future work and extensions of our system are presented in the conclusion.

2

The Pandore Environment

The PANDORE environment is shown in figure 1. It comprises a compiler, a machine-independent run-time and execution analysis tools including a profiler and a trace generator. The HPF front-end is built upon a Fortran 90 precompiler designed at GMD [START_REF] Jerid | HPF to C-Pandore Translator[END_REF]. The PANDORE run-time uses a generic message passing * E-mail: pandore(~irisa.fr

:: :'rt$ir s ~i i i ~! ~ i i!ii!ii!i:
library called POM (Parallel Observable Machine) [START_REF] Guidec | POM: a Virtual Parallel Machine Featuring Observation Mechanisms[END_REF]. This library offers limited but efficient services. It allows the same program to run on a wide range of distributed memory computers. The source language is a subset of High Performance Fortran. IIPF includes some "standard" extensions for data parallelism and data distribution into the Fortran 90 language. They permit to distribute arrays among virtual processor arrays. Examples of ItPF programs are shown in section 5, see [START_REF]High Performance Fortran Language Specification[END_REF] for a complete description of the language. When we designed the compiler, one of our objectives was to build a modular and extensible compiler in order to be able to integrate new techniques easily. The compiler has been written using the CAML language [2] that is a dialect of the functional language ML. The first prototype of the PhNDORE compiler relied on the well-known runlime resolulion technique [START_REF] Calahan | Compiling Programs for Distributed Memory Multiprocessors[END_REF] that has been proved correct [START_REF] Bareau | Correctness of Automated Distribution of Sequential Programs[END_REF]. This technique introduces masks and potential communications at the statement level. Consequently, it suffers from strong inefficiencies. So an optimized compilation scheme handling reductions and parallel loops with one statement has been added to the current version. A suitable run-time system completes these compilation techniques and permits to obtain efficient execution for a reasonably large class of programs.

Compiling Reductions and Parallel Loops. The compilation scheme is based on the decomposition of the arrays into blocks; it performs the restriction of iteration domains and the vectorization of messages. The mapping is taken into account at run-time. A complete description can be found in [START_REF] Le Fur | Andr& Static Domain Analysis for Compiling Commutative Loop Nests[END_REF]. In this section, we only describe our approach through the following (somewhat contrived) example where arrays A and B are partitioned into 8 blocks numbered from 0 to 7: For this parallel nested loop, the compiler produces a SPMD code that comprises a communication part followed by a computation part.

Communication Code Generation. First, the compiler performs a symbolic analysis on the iteration domMn and the array subscripts, allowing for the partitioning and the layout of the arrays, in order to construct the following system of attlne constraints:

{ 0<kA_<7 1<i< 1000 u=j
O<_kB<_7 i<__j<2,i+l v=i+j-2 500,kA<i< 500,~+499 500,kB<v <500,~+499 that defines the set of vectors (kA, kB, v, u, i, j) in which the array element B(u, v), located in the block number kB of B, is needed to perform the writings on the the block number kA of A. The above system of constraints defines a polyhedron P which is then projected along the i, j axis; this results in a new polyhedron Pij whose enumeration code can be computed by one of the methods described in [START_REF] Fur | Parcours de poly~dre paramdtrd avec l'dlimination de Fourier[END_REF][START_REF] Irigoin | Scanning Polyhedra with DO Loops[END_REF][START_REF] Collard | Construction of DO Loops from Systems of A]jine Constraints[END_REF]. This yields the nested loop (in pseudo-code): forM=O, 2 for kB = .~(0, 2 9 ~4 -1), mi~(5, 3 9 ~4 + 2) for v = max(500 * kB, 1000 * "/cA -2), rain(500 * kB + 499, 1500 * "kA + 1496) for u = max(div(v + 3, 2), v -998,-500 * kA -I-v -497),

.

~in(di~(2 9 ,, + 5, 3), -500 9 ~ + ,, + 2)

From this loop, the compiler generates the send code and the dual receive code for reference B(j, i + j -2). For the send code for instance, the compiler inserts send statements and appropriate masks in the loop so that the mapping of the blocks is taken into account at run-time: to owner_block(A, kA)

for kA = 0, 2 if myself ~ owner_block(A,
Although this loop contains masks, it is important to notice that these masks are evaluated at the block level and not at the iteration vector level as in the runtime resolution. Furthermore the (kA, kB)-loop do not scan the whole cartesian product 0..7 x 0..7 and the location of the first mask prevents from enumerating all the vectors described by the (kA, kB)-loop.

From the description of the elements B(u, v) to be sent, the run-time library routine RLR_send performs several communication optimizations. Direct communication is performed when possible: what is transferred in this case is a memory zone that is contiguous both on the sender and the receiver side, thus eliminating any need of coding/decoding or copying between message buffers and local memories. Message aggregation is also carried out and reduces the effect of latency by grouping small messages into a large message. Finally, the redundant communications, that may occur when several references to the same distributed array appears in the right hand side, are eliminated.

Computation Code Generation. The SPMD computation code is generated as follows: the compiler analyzes the iteration domain, the subscripts and the array partitioning for the reference A(i, j-i) to synthesize the set of constraints: 0 < kA< [START_REF]High Performance Fortran Language Specification[END_REF] 1 <i< 1000 kA<i<500..kA+499 i<j<2,i+l that defines the set of vectors (kA, i, j) where iteration vector (i, j) is such that reference A(i,j -i) writes in block ]cA of A. The enumeration code for the polyhedron associated with the previous system is then computed as in the communication code generation:

for kA = 0, 2 for i = max(500 * kA, 1), rain(500 9 "kA + 499, 1000)

forj=i, 2.i+1

From this nested loop, the computation code is finally generated by inserting an adequate mask so that the owner-writes rule is ensured at run-time:

for kA = O, 2 if myself = owner_block(A, M)
for i = max(500 9 "kA, 1), rain(500 * kA -4-499, 1000) forj=i, 2.i+1 A(i,j-i) := B(j,i + j-2)

As in the communication code, one can note that the number of tests performed by each processor is very small. First, the mask used to take into account the mapping at run-time is introduced at the block level and second, the outer M-loop does not scan the whole interval 0..7.

Management Scheme for Distributed Arrays

Representation of distributed arrays as well as accesses to elements of these arrays is a critical issue for overall performance of the produced code. In the PANDORE environment, arrays are managed by a software paging system. The run-time uses the addressing scheme of standard paging systems but is not a virtual shared memory: the compiler always generates communication when distant data are needed, so we do not need to handle page faults. The array management is based on the paging of arrays -not of memory:

tile multi-dimensional index space of each array is linearized and then broken into pages. Pages are used to store local blocks and distant data received. If data have to be shared by two processors, each processor stores a copy of the page (or a part of the page) in its local memory. Array elements are accessed through a table of pages allocated on each processor. The compilation technique ensures that accessed pages are up to date, hence the consistency between copies of array elements does not need to be handled at run-time.

Description. To access an element referred to by an index vector (i0,..., i~_ 1) in the source program, a page number and an offset (PG and OF) are computed from the index vector with the linearization function t; and the page size S: PG = C(io,...,i,~-l) div S, OF = Z;(i0,...,i,~-l) rood S. For a given distributed array, the parameters we tune for paging are the page size S and the linearization function L:. Time consuming operations are avoided in the computation of the tuple (PG, OF) but also in the evaluation of the flmction t; by introducing powers of two, turning integer division, modulo and multiplication into simple logical operations (shift, and mask). We first choose the dimension 6 in which the size of the blocks is the largest. Function s is the C linearization function applied to a permutation of the access vector that puts index number 5 in last position. The page size S is then defined by the following (se is the block size in dimension 5): if s~ is a power of two or dimension 5 is not distributed, S is the smaller power of two greater than s~; otherwise S is the largest power of two less than s~. Actually, an optimized computation of (PG, OF) is achieved by avoiding the explicit computation of the linear address Z;(i0,..., i,~-1): we express PG and OF directly as a function of the index vector, thus, when dimension ~ is not distributed, rood and div operations are removed. A more detailed description of this array management can be found in [START_REF] Mah6o | Distributed Array Management for HPF Compilers[END_REF].

Benefits. This management scheme leads to an efficient access mechanism. Access times remain very close to access times without index conversion and may be an order of magnitude faster than a "classical" index conversion involving a modulo or an integer division. The memory overhead induced by paging does not exceed a few percents for most distributions; it is almost entirely due to the tables of pages: when a page contains elements that have no equivalent in he original sequential space, or when just a part of a distant page is accessed in a loop, only a portion a the page is actually allocated. Apart from the performance aspects, paging distributed arrays offers several worthwhile characteristics. First, the scheme is always applicable. It is also independent of the analysis of the code: it only depends on distribution parameters, so no data re-arrangement or extra calculation is needed withing the scope of one distribution even if different loop and access patterns are concerned. The scheme is uniform: as far as accesses are concerned, no difference is made between local elements and distant elements previously received. Finally, the memory contiguity is preserved in the direction of the pages: contiguous elements of the original array are still contiguous in the local representation, facilitating direct communications and exploitation of caches and vector processors.

To our knowledge, management of distributed arrays have not been studied independently of compilation techniques. Several techniques such as the overlap, temporary buffers and hash tables are presently used in existing ItPF-like compilers [START_REF] Zima | SUPERB: A Tool for Semi-Automatic MIMD/SIMD ParMlelization[END_REF][START_REF] Tseng | An Optimizing Fortran D Compiler for MIMD Distributed-Memory Machines[END_REF][START_REF] Zima | Compiling for Distributed-Memory Systems[END_REF]. Although they have not been integrated in complete prototypes, other techniques that aims at packing local data to decrease memory overhead have been proposed [START_REF] Irigoin | A Linear Algebra Framework for Static HPF Code Distribution[END_REF][START_REF] Gilbert | Generating Local Adresses and Communication Sets for Data-Parallel Program[END_REF]. None of these methods gather all the characteristics discussed above.

5

Experimental Results

The compilation of several well-known kernels and larger applications have been tested with the PANDORE environment. Performance results are shown here for three kernels: Cholesky factorization, Matrix-matrix product and Jacobi relaxation; the description of the parallelization of a wave propagation application can be found in [START_REF] Andre | Parallelization of a Wave Propagation Appfication using a Data Parallel Compiler[END_REF]. The source code of the HPF subroutines are given. Apart from the distribution specification, they include minor modifications compared with the original sequential code. These modifications are to a large extent aimed at taking advantage of collective communications.

Measurements have been performed on a 32-node iPSC/2. The presented graphs show the speedup against the number processors for several input sizes (the indicated number is the value of N). Speedup is defined as the parallel time over the time of the original sequential program measured on one node. The obtained eificiencies are satisfactory, ranging from 85% to 95% on 8 processors and reaching around 80% on 32 processors for the largest data size. To handle very large applications, in particular applications that comprise multiple modules or necessitate intensive I/O, other techniques must be integrated to the existing environment. For this purpose, the joint study of redistribution, procedures and separate compilation is under way in the PANDOgE project.

Fig. 1 .

 1 Fig. 1. The PANDORE Environment

  kA) for ~B = max(O, 2 9 ~ -1), rain(5, 3 * kA + 2) if myself = owner_block(B, kB) for v = max(500 * kB, 1000 * kA -2), rain(500 * kB + 499, 1500 * kA + 1496) RZR_~ena {B(u, v) / u _> ma~(d~(v + 3, 2), ~ -998, -500 9 ~A + ~ -497) <_ mi,,(div(2 9 v + 5, 3),-500 9 ~ + v + 2))

  I-I,J) + B(I+~,J) + B(I,J-1)+ B0,J+~)) END DO END DO B(I:N-2, 1:N-2) = A(I:N-2, l:N-2) A,B,C) INTEGER, PARAMETER :: N = 512 REAL(KIND=8), DIMENSION(0:N-1,0:N-1) :: A, B, C !HPF$ PROCESSORS PROCS(4) !HPF$ DISTRIBUTE (BLOCK, *) ONTO PROCS :: A, C !HPF$ DISTRIBUTE (*, BLOCK) ONTO PROCS :: B INTEGER I,J,K REAL (KIND=8), DIMENSION (0:N-I) :: COLJ C(0:N-I,0:N-1) = 0.0 DO J=0, N-1 COLJ = B(0:N-I,J) DO I=0, N-1 DO K=0, N-i (2(i,J) .-= C(1,J) + A(I,K) * COLJ(INTEGER, PARAMETER :: N = 512 REAL(KIND=S), DIMENSION(0:N-1,0:N-1) :: A !HPF$ PROCESSORS PROCS(4) !HPF$ DISTRIBUTE (*, CYCLIC) ONTO PROCS :: A INTEGER I,J,K REAL (KIND=8), DIMENSION (0:N-I) :: COLK DO K=0, N-1 A(K,K) = SQRT(A(K,K)) DO J=Kq-1, N-1 A(J,K) = A(J,K) / A(K,K) END DO COLK(K+I:N-1) = A(KTI:N-1,K) DO J=K-t-1, N-1 DO I=J, N-1 h(I,J) = i(l,a) -COLK(I) * C, OLK(above described optimization techniques, the performances obtained on a series of numerical applications are already quite satisfactory even though enhancements can be made along several axis. For example, taking the mapping of the blocks into account at compile time will allow us to suppress masks in the loops generated in the communication and the computation codes. Moreover, we plan to enlarge the subset of HPF compiled, especially by adding alignments and nested subroutines calls.

We think that these improvements will contribute to a better maturity for data parallel compilers and so automatic code generation for distributed memory parallel architectures will become a realistic means of programming these architectures for application users.


Yves Mah

Jean-Louis Pazat

Distributed Array Management S c heme for Data-Parallel Compilers

Keywords: Distributed memory parallel computers, compilation, data-parallel languages, HPF, runtime, paging, memory management

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In order to alleviate the task of programming Distributed Memory Parallel Computers, new features have been added to sequential programming languages such a s F ortran. In the eld of scienti c programming, two main axes are currently followed. The rst one uses explicit parallel constructs and loop partitioning, it relies on a shared virtual memory 2] the second one is based on a userde ned data distributions which are used as a guideline to generate communicating processes 7].

In recent y ears, many projects have focused on the data distribution approach and it has been demonstrated that \aggressive" optimizing compilers and e cient r u n time systems are mandatory to achieve reasonable speedups. Most compilers allow the user to specify a decomposition of arrays and use the owner write rule 4] to distribute the code. This distribution can be done using the runtime resolution technique in which e a c h statement of the source program is guarded and where communication is performed elementwise. This scheme is always applicable but rather ine cient, so that many compilers integrate optimization techniques for compiling loops. Roughly speaking, 1 these techniques aim at reducing iteration domains and performing vectorized communications 14, 12, 1 1 , 1 0 ] .

However, runtime resolution as well as optimization techniques require a speci c and e cient runtime system to allocate local parts of arrays and to perform e cient c o m m unications.

In this paper we present a new method for e ciently managing distributed arrays (allocation and access) in parallelizing compilers based on data distribution. This paper focuses on the local management of blocks and temporary storage for distant v alues. Communication optimizations such as message coalescing and vectorization are supported by this array management but are not addressed here.

The paper is organized as follows: next section discusses the essential requirements for managing distributed arrays. Section 3 details the page-driven array management proposed. Section 4 presents in more details the implementation of this array management whose overall performance is presented in section 5. Future work is discussed in the conclusion.

2 Management of Distributed Arrays

Key issues

The aim of a distributed array management is to de ne the way distributed arrays are stored in the local memories and the way e l e m e n ts of these arrays are accessed. These tasks involve both the compiler and the runtime support. It is clear that a trade-o must be achieved between the speed of accesses and the memory overhead induced from the array representation. The extreme solution consisting in allocating the entire array o n e a c h processor is obviously not applicable. Conversely, minimal allocation (typically achieved by translating an array declaration A(N) i n to a local allocation A(N=P) w h e r e P is the number of processors) associated with global-to-local index conversion involving several costly operations such a s mod and div must be avoided.

In addition to the \classical" concerns that are the memory use and the speed of accesses, we de ne in the following some properties that we claim to be useful for a distributed array management.

Uniformity

Two kinds of data are to be considered on a given processor: local elements { i.e. elements assigned to it by the distribution { and received elements that are temporarily stored in the local memory after communication with another processor. We s a y that an array management s c heme is uniform when the representation as well as the access mechanism associated with a distributed array m a k es no distinction between local and received elements. The main advantage of a uniform scheme is that the compiler does not have to separate purely local computation (involving only local data) from computation that needs distant data. Indeed, even if such a separation is sometimes possible at compile-time, it may induce an important code fragmentation in the case of multiple right-hand-side references performing the separation at runtime brings about costly ownership computation.

Furthermore, de ning a non-uniform scheme that gives too great an importance to local elements, at the expense of received elements is not a good solution because it is likely to be harmful not only for accesses but also for communication of non-local data.

Independence

An independent array management is only de ned from distribution parameters and in particular does not depend on the code itself nor on compiler analysis on the code. A distributed array management s c heme independent from the compilation scheme facilitates the coexistence of di erent compilation techniques. On the contrary, i f t h e c hoice of a representation is guided by the analysis of a program part (typically a loop), it may happen that several layouts (and associated access methods) are used within the scope of one distribution, possibly necessitating data rearrangement or additional computation at runtime. Moreover, this independence facilitates the use of di erent compilation techniques within a code fragment that contains several loops. One way t o a c hieve this independence (as well as uniformity) is to consider that only global indices appear in the generated code.

Contiguity Preservation

Another useful property concerning the layout of distributed arrays is the conservation of memory contiguity. Indeed, if contiguous elements of the original array are still contiguous in the local representation, it makes it possible to take a d v antage of direct communications (in this case no copying nor packing/unpacking is needed between local representations and communication bu ers), vector processors, target code optimization and better cache behavior.

Related Work

To our knowledge, management of distributed arrays have not been studied independently from compilation techniques in existing HPF compilers.

The rst technique of storage for distributed arrays, the overlap 13] has been implemented both in the Vienna Fortran Compilation System 14] and in the Fortran D compiler 8, 12]: a single subarray is allocated for local data as well as for received data. This technique provides uniform accesses and preserves memory contiguity but it can be applied to a restricted number of distributions and access patterns and may lead prohibitive allocation when distant data location is not close to the local partition. In this case, the Fortran D compiler may select an alternative storage method (bu ers) for received values if it can separate purely local computation and computation needing received values. In Vienna Fortran and in an extended version of Fortran D, a speci c management related to loops with irregular array accesses is performed through the use of the inspector/executor technique: the local partition is dynamically extended during the inspector phase so that uniform indirect accesses can be used during the executor phase 6, 3 , 1 4 ].

Other array management s c hemes, closely related to compilation techniques, have been proposed. Their common characteristic is that they try to minimize memory overhead. Among them, in the compilation scheme de ned by Ancourt et al 9], local elements and temporaries are packed according to the array distribution and alignment, the loop bounds and the array subscripts by changing the basis of the original index space. Accesses to elements are performed in a non-uniform way with index conversion evaluating a ne functions and possibly integer division. Chatterjee et al. 5] propose an access mechanism for local elements based on a Finite State Machine (FSM). These elements are accessed by executing a FSM that has to be computed at runtime for each l o o p nest even if the same distribution applies.

All of these methods not only take i n to account the array distribution parameters but necessitate also a static analysis of code fragments (loop bounds and array subscripts) in order to de ne the layouts of the local arrays and the associated access mechanisms. Therefore, independence from the compilation scheme is not achieved by these systems.

A P age-driven Array Management

We p r e s e n t here a new management s c heme for distributed arrays based on software paging. This management is designed in order to achieve e cient accesses while avoiding unacceptable memory overhead. It also aims at satisfying the properties of uniformity, independence and contiguity preservation aforementioned. In the following, we will consider only direct HPF distributions, the mechanisms described can be easily extended to aligned distributions (as a rst step, by applying paging to templates).

Principle

The page-driven data management w e propose follows the main addressing scheme of classic paging systems for memory management. In such systems, logical memory space is broken into groups of contiguous elements (pages). Pages have a xed predetermined size. A hardware support divides a logical address in two parts: a page number and a page o set. The page number is used as an index into a page table that contains the base address of each p a g e i n p h ysical memory. This base address is combined with the page o set to de ne the physical memory address. If the page size is S, a logical address produces a page number PG and an o set OF by PG = div S and OF = mod S. If the logical address space is larger than the physical address space, virtual memory management features may be added. In this case, accessed pages may not be present all the time in physical memory but temporarily loaded from a secondary storage system by a swapping device.

As for our concern, we manage variables |i.e. distributed arrays| and not memory our aim is not to build a shared virtual memory. Moreover, we only consider compiler-generated code, hence we s t a y at the software level rather than relying on hardware components. Contrarily to system-level paging, the notion of page fault is here irrelevant because all distant accesses are solved by prior communications. Besides, data are not necessarily communicated page-wise.

The original address space is multidimensional therefore we a p p l y a m ulti to one-dimensional transformation before splitting the resulting space into pages. This allows us to de ne a speci c access mechanism for each distributed array, in particular the page size may be di erent f o r e a c h array.

Paging Distributed Arrays

We de ne a representation and its access mechanism for each distributed array b y a couple (L S ).

The multidimensional index space of a given array is linearized by a function L. The linear address space obtained is split into pages of xed size S. A processor stores only those pages that contain at least one element a s s i g n e d t o i t b y the distribution or one received element. Depending on the distribution of the array, L and S, a page may be possessed by one or several processors.

One of the main advantages of this method is that accesses to local and received elements are performed the same way. Indeed, as far as accesses are concerned, a processor acts as if the entire array w as allocated locally, no matter if the element it needs to access is truly local or has been received from another processor. The di erence between pages containing local elements and pages containing received elements lies in the way they are allocated and lled, not in the way they are accessed. A tuple (PG,OF ) is computed from the initial index vector (i 0 : : : i n;1 ) w i t h t h e linearization function L and the page size S: PG = L(i 0 : : : i n;1 ) div S OF = L(i 0 : : : i n;1 ) mod S

A table of pages is stored on each processor. It indicates the base address of each page present i n local memory. The o set is added to this base address to obtain the exact location of the element. The page partitioning is also used for computing owners of elements. A similar table, present i n the local memory of each processor, stores for each p a g e t h e n umbers of the physical processors that own this page.

Tuning Parameters

For a given distributed array, the parameters we can tune for paging are the page size S and the linearization function L. T h e v alue of these parameters are de ned in order to achieve good performance in terms of time and memory space.

As speed of access is our prior motivation, time consuming operations (division, modulo and multiplication) are avoided in the computation of the tuple (PG,OF) but also in the application of function L. This is achieved by i n troducing powers of two, turning integer division, modulo and multiplication into simple logical operations. Moreover, the array decomposition can be taken into account when xing the actual value of S and L. I n tuitively, w e c hoose S and L so that pages \follow" the blocks, and are owned by as few processors as possible.

For a more formal de nition, let us consider the following HPF array distribution: REAL V(0 : h 0 ; 1 0 : h n;1 ; 1) !HPF$ DISTRIBUTE V(CYCLIC(s 0 ),: : : ,CYCLIC(s n;1 )) and the access to an element o f V noted V(i 0 i n;1 ). This distribution decomposes the array V into rectangular blocks of size s 0 s n;1 . W e consider the most general distribution directive CYCLIC(k). Note that BLOCK(k) and CYCLIC(k) distributions are strictly equivalent a s f a r a s decomposition is concerned.

Prior to the de nition of S and L, w e c hoose a particular dimension , the dimension in which the block size is the largest. If there are several such dimensions, the one corresponding to a nondistributed array dimension or a block size equal to a power of two i s c hosen. The page size S is then given by: if s = h or s = 2 then S = sup (s ) else S = inf (s ) where sup (x) (resp. inf (x)) extends an integer to the smallest (resp. largest) power of two greater (resp. less) than or equal to x.

L is the C linearization function for multidimensional arrays applied to a permutation of the index vector. This permutation puts the index corresponding to dimension in last position.

Moreover, the array dimensions (coe cients of L) are extended to the next power of two. L is de ned by L(i 0 : :

: i n;1 ) = n;1 X k=0 i 0 k n;1 Y l=k+1 h 0 l !
where i 0 k is the k th access index after permutation, i.e: i 0 n;1 = i 8k 2 0 : : : ;1 i 0 k = i k 8k 2 : : : n ;2 i 0 k = i k+1 and h 0 k is the extended size of the array i n t h e k th dimension, i.e:

h 0 n;1 = sup h S S if n > 1 h 0 0 = h 0 if > 0 else h 1
8k 2 1 : : : ;1 h 0 k = sup (h k ) 8k 2 : : : n ;2 h 0 k = sup (h k+1 )

Here are two examples of de nition of S and L rst when there is one non-distributed dimension and second when all the dimensions are distributed:

REAL A(0:199, 0:99, 0:50) !HPF$ DISTRIBUTE A(CYCLIC(5), , CYCLIC(10)) ) S = 1 2 8 L(i j k) = ( 6 4 128)i + 1 2 8 k + j REAL B(0:499, 0:199) !HPF$ DISTRIBUTE B(CYCLIC(100), CYCLIC(10)) ) S = 6 4 L(i j) = 5 1 2 j + i

Optimizing the Computation of (PG,OF)

Unlike with a classic paging mechanism, the explicit computation of the linear address L(i 0 : : : i n;1 ) before its splitting into (PG,OF) is not mandatory because we do not rely on a hardware support that needs a memory address. Besides, this intermediate result may lead to unnecessary operations as in the following example:

REAL A(0:99, 0:199) !HPF$ DISTRIBUTE A(CYCLIC(10), ) ) S = 256 L(i j) = 2 5 6 i + j

The page number and the o set will be obtained by PG = ( 2 5 6 i + j) div 256 OF = (256i + j) mod 256 These expressions could obviously be simpli ed in PG = i and OF = j. T o m a k e the simpli cations clearly visible, we express directly PG and OF as a function of the index vector.

page(i 0 : : : i n;1 ) = ( PG OF) with PG = n;2 X k=0 i 0 k n;1 Y l=k+1 np 0 l ! + i 0 n;1 div S OF = i 0 n;1 mod S where np 0 k is the number of pages in the k th dimension after permutation: np 0 n;1 = h 0 n;1 S 8k 2 0 : : : n ;2 np 0 k = h 0 k When dimension is not distributed, that is to say w h e n h = s , index i 0 n;1 (i.e i ) i s a l w ays less than or equal to S, div and mod can be removed:

Page Ownership

Each processor stores a table of owners that indicates, for each page, the number of the physical processor that owns this page. This table can be lled using the function owner(PG,OF ) that returns the owner of an element. owner(PG OF) = map page ;1 (PG OF) Function page ;1 , the reverse function of page, returns the index vector corresponding to a page number and an o set.

page ;1 (PG OF) = ( i 0 : : : i n;1 ) with i = S ; PG mod np 0 n;1 + OF 8k 2 0 : : : ;1 i k = i 0 k 8k 2 +1 : : : n ;1 i k = i 0 k;1 8k 2 0 : : : n ;

2 i 0 k = PG mod n;1 Y l=k np 0 l ! div n;1 Y l=k+1 np 0 l !!
Function map(i 0 : : : i n;1 ) associates a physical processor number with an index vector. It can be easily computed for each mapping that can be expressed in HPF. We d o n o t g i v e here the general formula as it cannot be induced from the HPF norm and depends on implementation choices. As an example if the abstract processor array is of size (p 0 p n;1 ) and the numb e r o f p h ysical processors is P, the mapping function may be the following:

map(i 0 : : : i n;1 ) = n;1 X k=0 (i k div s k ) n;1 Y l=k+1 p l !! mod P
In the case data elements are replicated (by application of alignment directives), this function could return a set of processors. The de nitions adopted for S and L allow t h e n umb e r o f o wners of a page to be less than or equal to two. If the owner of a page is always unique, any v alid value of OF can be used for determining the owner of a page. In the case the owner of a page is not unique, we can compute OF lm , the o set from which t h e o wner changes. In this case, the table of owners stores for each page, the two processor numbers plus the limit OF lm :

OF lm = i f ' < S then ' else 0 with ' = ( ( PG mod np 0 n;1 ) (s ; S)) mod s 4 Implementation A full implementation of the data management m e c hanisms described above has been realized within the Pandore environment 1 ]. In the Pandore language, the scope of array distributions is con ned within procedures referred to as distributed phases.

Management of tables, pages and accesses to array elements are shared out among the compiler and the runtime library. As all the tables and pages are needed only during the execution of a distributed phase (no inter-phase analysis is performed at this time), the entire memory space allocated is freed at the end of the phase.

Tables and Pages

All the information needed to ll the tables of owners and the tables of o set-limits is known at compile-time these tables could therefore be statically de ned. However, in order not to lengthen the size of the generated code, the compiler produces functions that allocate and ll the tables at runtime, at the beginning of each distributed phase. For each distributed array V , a table of owners TO V i s d e n e d . I f a p a g e m a y be possessed by t wo processors, three tables are needed: the table of the owners of the rst part of pages TO1 V, the table of the owners of the second part of pages TO2 V and the table containing the o set-limits TL V.

The runtime library is also in charge of allocating and lling the tables of pages and pages themselves. The tables of pages and pages that contain local elements are allocated at the beginning of the distributed phase. The management of pages containing received elements depends on the compilation scheme. Basic operations provided by the runtime library are the page allocation and the placement of elements (single elements or segments) into pages.

Accesses

It is clear that the part of the access process that is done at compile-time must be as large as possible. The compiler translates a reference to an array element V I], where I is an index vector, into a call to a runtime macro access(desc V, PG, OF) where PG and OF are expressions of I. All constant subexpressions have been computed and the optimization described in section 3.4 has been performed. As expected, these expressions contain only additions and constant logical shifts and maskings. The work that remains at runtime is therefore to evaluate the expressions and use the table of pages associated with V (TP V) to produce the right reference. This can be noted by the C expression *(TP V PG]+OF). The runtime library contains cpp macros that prevent from the computation of the address of the page table corresponding to V , s o w e can actually generate this code.

Owner

Determining the owner of an element V I] is carried out a similar way. The compiler generates a call to a runtime macro owner(desc V, PG, OF). An access to a table TO V PG] is su cient a t runtime to nd the processor number in the case the owner of a page is unique. If a page may b e possessed by t wo processors, a call to a slightly di erent macro is produced. The execution of this macro will issue a comparison between OF and the o set-limit corresponding to page It is quite obvious that the executed code for distributed accesses involves only few basic operations that generate a very small overhead and may e v en be more e cient thanks to better optimizations.

In the experiment whose results are reported in table 1, we measured the time taken by several kinds of read accesses: t s : a reference to an element a s i t m a y appear in a sequential program t p : a call to the macro that uses the paged access mechanism t b : a call to a macro that uses a block-oriented access mechanism1 . The array i s a t wo-dimensional array of oats reported times are in s. Best and worst cases have been considered, depending on whether the sizes of the array w ere powers of two or not. Experiments have been carried out on a SparcStation 2, on a node of the iPSC/2 and on a node of the Intel Paragon XP/S.

Likewise, the determination of the owner of an array e l e m e n t requires only a few simple operations, so its cost remains very low. It is also preferable to exploit the page decomposition, although it seems to be more natural to base the computation of the owner of an element on the computation of the corresponding block n umber.

The price to pay for speed of access and speed of ownership computation is the need for a larger amount of memory. O v erhead is only due to tables because no additional space is required for pages. When a page contains elements that will never be accessed because of the extension of array dimensions, or because the page is shared by t wo processors, only the potentially accessed part of the page is actually allocated. A translation of the corresponding pointer in the table of pages is performed if the end of the page is allocated.

The memory overhead due to tables is directly linked to the number of pages, which is in general at least of an order of magnitude less than the size of the array. T able 2 gives memory requirements for a few common distributions of arrays on 32 processors. For each distribution, we indicate the total number of pages, the theoretical minimal memory space required on each processor, the actual space allocated for tables on each processor and nally the overhead as compared with the minimal partition. Memory needs are expressed in bytes. It can be noticed that replacing some

Integration in the Pandore Environment

The page-driven management for distributed arrays has been integrated in the Pandore environment and is used with the two compilation schemes of the compiler. The basic compilation scheme, that relies on a runtime resolution technique, can be applied to every input program. The optimized scheme is based on integer programming and linear algebra results it performs an analysis of parallel loops 10]. We p r e s e n t in table 3 the results of the execution of a Red-Black Successive O v er-Relaxation algorithm run on a 1024x1024 matrix of oats. Times have been measured on the Intel iPSC/2 for the two compilation schemes. A comparison is made between a block-oriented array management and the page-driven management. The table shows the speedup obtained on P processors for each pair (compilation scheme, array management). The use of the page-driven management clearly improves performances of codes generated cording both compilation schemes. The joint use of the optimized scheme and the page-driven array management leads to satisfactory performances (e ciency of around 75% for 32 processors) in spite of the unfavorably high ratio of memory operations to computation of the Red-Black algorithm.

Conclusion

Management of distributed arrays is a crucial point for obtaining good performances when using data parallel compilers. For this purpose, we h a ve proposed a new scheme based on parameterized software paging that proved e cient in an existing compiler. This management handles local and received data in an uniform way and it is independent from the optimization techniques used in compilers. Moreover, it avoids using multiple representations of the same array in di erent parts of a program and maintains some regularity i n l o c a l l a youts. The page-driven array management also seems to be appropriate for irregular computations and could be used together with the inspector/executor technique.

We are currently comparing our management s c heme with shared virtual memory systems. One of our objectives is to nd out which features should be added to existing shared virtual memory systems so they can be e ciently targeted by data-parallel compilers.

  is the result of these optimizations for the two examples presented in the previous section: REAL A(0:199, 0:99, 0:49) !HPF$ DISTRIBUTE A(CYCLIC(5), , CYCLIC(10)) )

Table 2 :

 2 Memory overhead for some common distributions

block sizes (or array dimensions) by p o wers of two notably decreases the memory overhead. We believe t h a t o verall memory requirements remain acceptable when considering most commonly used distributions.

Table 3 :

 3 Comparison between block-oriented and page-driven managements

	Basic scheme Optimized Scheme P Block Page Block Page
	4 0.29 0.68 0.86 8 0.32 0.77 1.39 16 0.36 0.82 2.14 32 0.37 0.85 3.70	3.84 7.18 12.78 23.72

This mechanism was used in a previous version of Pandore it performs at runtime a modulo and an integer division to nd the block n umber and the o set in the block.


Olivier Delestre

François James

García De Galdeano

SIMULATION OF RAINFALL EVENTS AND OVERLAND FLOW

Keywords: shallow water equations, finite volume schemes, well-balanced schemes, hydrostatic reconstruction, friction laws, rainfall hydrograph, analytical solution, dry/wet transitions AMS classification: 76M12, 74G05, 74G75, 35L65, 20C20

We are interested in simulating overland flow on agricultural fields during rainfall events. The model considered is the shallow water system (or Saint-Venant equations) without infiltration, complemented with a friction term. In this context, we definitely have to cope with dry/wet interfaces and water inflow on dry soil. We present a simplified one-dimensional model, discretized with a well-balanced finite volume method, and we describe the specific additional features needed to deal with dry/wet transitions and steady-state solutions due to topography and friction. The method as well as the choice of the friction term are tested and discussed both on analytical solutions and experimental results.

Introduction

Rain on agricultural fields can yield to overland flow. This flow may have some undesirable effects. At the field scale, we can have soil erosion and pollutant transport. Downstream the fields, roads and houses may be damaged. To prevent these effects, control measures can be taken, such as grass strips. But one must know how the water is flowing in order to place efficiently these developments. In the spirit of [START_REF] E | Overland flow and infiltration modelling for small plots during unsteady rain : numerical results versus observed values[END_REF][START_REF] F | A numerical method for simulating discontinuous shallow flow over an infiltrating surface[END_REF], we try to model these phenomenon by using the shallow water (or Saint-Venant) equations. Efficient numerical simulations are of great help in this context, because field measurements, such as velocities or water heights, are very difficult to obtain, especially during the rain event, which is quite unpredictible.

The aim of this paper is not to give a complete account on the problem, which has to be thaught of as a multi-scale problem: one has to deal with roughness induced at the decimeter scale (e.g. by furrows on agricultural surfaces), flows at the scale of ten square meters, which is the scale of the numerical topography data, and also the agricultural field itself, whose surface is of the order of the hectare. We give here a short review of the shallow water equations, with emphasis on some specific aspects in this context. Namely, since the rain is an intermittent phenomenon, we definitely have to cope with dry/wet transitions, a problem analogous to the vacuum apparition in gas dynamics. More classically in shallow water problems, we have to take into account carefully the interactions between the soil topography and the friction of water on the soil, which eventually lead to steady-state solutions that have to be computed accurately.

For this introduction to the topic, we deliberately use a simplified model, firstly by considering one-dimensional flows. This is enough to understand the ideas of the numerical methods, which can be developed in two space dimensions on a rectangular mesh. Next, from a more practical viewpoint, we neglect importants phenomena, which deserve a complete modelling: infiltration and soil erosion. Infiltration appears as a supplementary source term in the shallow water equations, and can be treated quite easily, when a relevant model is chosen. Erosion is a much more complex problem.

We begin by a short review of the shallow water system, recalling a few basic properties. Next, we describe numerical methods adapted to the situation, in particular we discuss briefly the discretization of the friction terms. Finally, we give several illustrations of the results. First we justify the choice of the method by comparison with analytical solutions. Next, we show an attempt of recovering experimental results, with a clear evidence that the choice of the friction laws is not obvious. The last section is devoted to an unstability phenomenon wich occurs when perturbating steady-state solutions (with rain for instance): the so-called roll-waves. §1. Model

The model we consider here are the so-called shallow-water equations, which are convenient for small heights of water, according to the following scheme

x O h(t,x) z(x) u(t,x) z+h z
The unknowns are here the velocity of the water u(t, x), and its height h(t, x). The shape of the bottom is also called the topography, it is a given function z. For our specific application, the model has to be complemented by taking int account friction on the soil and rain. Therefore the equations are

∂ t h + ∂ x (hu) = R(t), ∂ t (hu) + ∂ x hu 2 + gh 2 2 = -gh ∂ x z + S f , (1)
where g is the gravity constant, R(t) the rain intensity, assumed constant in space, and S f (h, u) the friction term. Notice that infiltration in the soil can be accounted by a source term in the first equation like R(t) -I(t, x), where I is a given function. We shall denote by q = hu the water flow, or discharge. The typical practical configuration we consider is a channel with finite length L, so that the system must be set on the interval ]0, L[, and complemented with boundary conditions at inflow and outflow we do not detail here, see an example in Section 3.

Concerning the friction term, it is a given function of h and u, two examples widely used in hydrology (see for instance [START_REF] E | Overland flow and infiltration modelling for small plots during unsteady rain : numerical results versus observed values[END_REF][START_REF] F | A numerical method for simulating discontinuous shallow flow over an infiltrating surface[END_REF][START_REF] Md | Steady open channel test problems with analytic solutions[END_REF][START_REF] Md | Analytic benchmark solutions for open-channel flows[END_REF]) are the Manning and the Darcy-Weisbach friction laws, which are given respectively by

S f = - k 2 u|u| h 4/3 = - k 2 q|q| h 10/3 , S f = - ku|u| 8gh = - kq|q| 8gh 3 , (2)
where k > 0 stands for the roughness coefficient. Both laws are derived from empirical considerations, in particular in the context of pipelines. The problem of their relevance in the present context of overland flow is difficult.

The system can be rewritten in a more compact form by setting

U = h q , F(U) = q q 2 /h + gh 2 /2 , B = R -gh ∂ x z + S f .
We obtain therefore

∂ t U + ∂ x F(U) = ∂ t U + F ′ (U)∂ x U = B.
The system is by definition hyperbolic if the matrix F ′ (U) admits a basis of eigenvectors with real eigenvalues, strictly hyperbolic if the eigenvalues are distinct. An easy computation shows that the shallow water system is strictly hyperbolic provided h > 0, with eigenvalues λ -(U) = ugh, λ + (U) = u + gh. When h = 0, the system is no longer hyperbolic, actually it is rather meaningless, since h = 0 means that there is no water, so that the velocity u cannot be defined. This is exactly the problem of the vacuum in the Euler equations of fluid mechanics, and leads to severe numerical problems, which cannot be avoided in our context since we consider rain on dry soils. At this point, we introduce an important quantity, the so-called Froude number

Fr = u gh . (3)
This dimensionless number plays the same role as the Mach number in fluid mechanics, and allows to classify the flows: -Fr < 1 subcritical flow, as in a river (corresponding to subsonic flow in fluid mechanics); -Fr > 1 supercritical flow, as in a torrent (subsonic flow); -Fr = 1 critical flow (transonic flow).

The differences between these flows can be easily experimented by observing the surface waves obtained by throwing a stone in a river. §2. Numerical method

The shallow water system is discretized by a finite volume method on a fixed time-space grid. A time step ∆t > 0 and a space step ∆x > 0 are fixed, we set x i = i∆x, and the interval ]x i -∆x/2, x i + ∆x/2[ will be referred to as the cell i. The finite volume scheme can be written in a compact form as

d dt U i + 1 ∆t (F i+1/2 -F i-1/2 ) = S i , (4)
where the vector U i is an approximation of the conservative variables in the cell i, F i+1/2 is the numerical flux at the interface between cells i and i + 1, and S i a discretization of the source term. Boundary conditions are treated by the method of characteristics (see [START_REF] B | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF]). The scheme is completely determined once the numerical flux and the source term discretization have been fixed. These choices are not independant one from the other. Indeed it is well-known that source terms in hyperbolic systems of conservation laws give rise to serious problems. The main difficulty is to find schemes that preserve equilibria (steady-states solutions). In system (1), the main problems are due to -topography: pools, lakes; -friction terms: balance between kinematics and friction. The rain source term can be treated by a second-order accurate Strang type splitting.

Schemes that preserve equilibria are known as well-balanced schemes. The strategy to obtain such schemes consists in choosing first a consistent numerical flux for the system without source terms. Next, a correction is given to take into account equilibria. The reader can find all the details and a large bibliography in the book [START_REF] B | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]. We merely give a sketch of the method here, with emphasis on the problem of friction. The numerical flux is the so-called HLL flux, and the order 2 is obtained in space by a MUSCL type reconstruction, in time by Runge Kutta (Heun) (see [START_REF] B | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] for details). Notice that dry/wet transitions imply a specific reconstruction for the water height, not only for the velocity as usual (see [START_REF] A | Modélisation hyperbolique et analyse numérique pour les écoulements en eaux peu profondes[END_REF]).

First we consider the equilibria for topography. They are given by

hu = Cst, u 2 /2 + g(h + z) = Cst.
However a complete resolution of these equations would lead to a far too time consuming scheme. Thus, following [3, 1, 2], we limit ourselves to the equibria at rest:

u = 0, g(h + z) = Cst.
This procedure is known as the (second order) hydrostatic reconstruction, and it turns out to give good results at an acceptable numerical cost. We refer the reader interested into details to the preceding references. Now we turn to friction terms, which can be treated by two different means. The first one aims at building a well-balanced scheme for friction as well as topography, is the apparent topography method, introduced by [START_REF] B | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]. It consists in building an modified topography z app which takes into account the friction, as follows:

z app = z -b, with ∂ x b = S f .
We proceed then exactly as before, with this new topography (detailed computations for the friction laws (2) can be found in [START_REF] D | Simulation of Rain-Water Overland-Flow[END_REF]). This gives rise to a scheme which computes neatly equilibrium states, but is not completely satisfactory on transition solutions, as we shall see in the next section.

Therefore we turned to a splitting method, and we chose the semi-implicit treatment proposed in [START_REF] B | Boundary conditions for the shallow water equations solved by kinetic schemes[END_REF], not only because it preserves steady states at rest, but also for its stability. For the Darcy-Weisbach friction law (2)-right, it writes

q n+1 i + f |q n i |q n+1 i 8h n i h n+1 i ∆t = q n i + ∆t ∆x i (F i+1/2G -F i-1/2D ),
where the right-hand side is nothing more than the discharge obtained at each step of the second order in time Runge-Kutta reconstruction. Notice also the simplicity of the method, which gives an explicit value for q n+1 i . Now we illustrate these ideas on a set of analytical solutions. §3. Analytical solutions Here we present briefly an adaptation to the 1-d case and our friction laws of an idea presented in [START_REF] Md | Steady open channel test problems with analytic solutions[END_REF][START_REF] Md | Analytic benchmark solutions for open-channel flows[END_REF] for pseudo two dimensional cases. At steady states, we have ∂ t h = ∂ t u = ∂ t q = 0, thus the mass-conservation equation gives q = cst and we get the equation

∂ x z = q 2 gh 3 -1 ∂ x h + S f (q, h) (5)
where S f (q, h) depends on the friction law chosen, for instance [START_REF] A | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. For any given value of the constants k and q, once we are given an explicit expression for h(x), then formula (5) allows us to compute the topography corresponding to this steady state and this water height.

Other friction laws can of course be chosen.

As an example, we consider a channel of length 1000 m, with a specified water height h(x) given by

h(x) = 4 g 1/3        1 + 1 2 exp        -16 x 1000 - 1 2 2               .
The friction model is the Manning law, with roughness coefficient k = 0.033. The topography is calculated iteratively thanks to [START_REF] D | Simulation of Rain-Water Overland-Flow[END_REF]. To make use of the shallow water system, we have now to impose boundary conditions. Since the flow is subcritical both at inflow x = 0 and outflow x = 1000, we have to impose the value of one quantity at inflow and one at outflow. We choose to put a discharge of q = 2 m 2 /s at inflow and a water height corresponding to the value of h(1000) downstream. We first compare the results obtained by the apparent topography and the semi-implicit scheme in preserving the equilibrium state. It turns out that both methods preserve correctly the steady state along time, as is evidenced by fig. 1.

Since for our application we are particularly interested in non-stationary solutions, we have considered an initially dry soil and the upstream discharge q = 2 m 2 /s, and computed the unsteady solution up to equilibrium. Both methods (apparent topography and semi-implicit treatment) converge towards the steady state, with slightly better results with the apparent topography method. However, before the steady state is reached, we have a wet/dry transition (fig. 2). We note that the apparent topography method is not adapted to this transition: we In figure 3, two more examples of computation of steady states are displayed, both with sub-and supercritical inflows and outflows, and using the semi-implicit method. The numerical scheme deals in particular with transition from one regime to the other, including hydraulic jumps (fig. 3-right). For this test, dry/wet transitions are involved, since on the one hand there is no rain on the last 5 cm of the flume, on the other hand rain falls on a dry soil. The measured output is an hydrograph, that is a plot of the discharge versus time (see fig. 5). The mathematical model for this ideal overland flow is the following. We consider a uniform plane catchment whose overall length in the direction of flow is L. The surface roughness and slope are assumed to be constant in space and time. The friction law is the Darcy-Weisbach one. We consider a constant rainfall excess such that

R(x, t) = I for 0 ≤ t ≤ t d , 0 ≤ x ≤ L, 0 otherwise,
where I is the rainfall intensity and t d is the duration of the rainfall excess. First we compute some explicit "naive" analytical solution to the problem. We notice that three phases can clearly be identified on the hydrograph: a first non-steady step at the beginning of the rainfall event, then a steady-state and lastly another non-steady step when rain stops. The first and the second step solutions can be computed explicitly, and the "naive" solution is obtained by assuming a simple concatenation of the two parts (we refer to [START_REF] D | Simulation of Rain-Water Overland-Flow[END_REF] for the detailed computations). At first we compare numerical results with the analytical "naive" solution. Once again, with (fig. 6-a) we show that with the apparent topography method, we get a peak on the discharge downstream that we do not get far from this transition. With the semi-implicit method, we do not have this peak (fig. 6-b). This treatment gives good results close to the "naive" exact solution. The hydrograph is well calculated (fig. 6-b), notice here the computed hydrograph at the middle of the flume, a quantity hardly accessible by experiment. Next, we propose a comparison between experimental measurements and numerical simulation (fig. 5), obtained with the Darcy-Weisbach friction law. We obtain a reasonable agreement, but it turns out that it is impossible to fit correctly the shape of both the increasing and decreasing parts of the hydrograph. This indicates clearly that the model has to be modified, for instance by choosing alternative friction laws, but this is beyond the scope of this paper. §5. Roll waves This section is devoted to some examples of the so-called "roll-waves", a phenomenon which results from the competition between topography and friction. Several steady regimes turn out to be unstable, a slight perturbation generating a periodic travelling wave with shocks (hydraulic jumps). In ref [START_REF] Q | The numerical study of roll-waves in inclined open channels and solitary wave run-up[END_REF], Que and Xu gather a set of explicit computations in the simple case of a constant steady states in inclined open channels with constant slope. They provide a precise analysis for the linear stability, proving in particular the following criterion: the initial constant state is linearly stable if and only if the Froude number (3) is smaller than 2.

We recover here these results, using the semi-implicit scheme described above, together with hydrostatic reconstruction. The initial height of water is different for each case, but the amplitude of the perturbation is the same. The "final states" showed here are computed at time t = 200s., since it turns out that the solution is stabilized at this time. All cases are perfectly computed, the convergence rates for different values of the Froude number are given in figure 8.

Comparisons between the initial perturbation and the final state are displayed in fig. 7. For Fr = 2, the initial state is supposed to be exactly stable, the smaller amplitude of the final result is due to the numerical diffusion. Notice the nonlinear effects (fig. 7, top right). For Fr < 2 (top left), the initial perturbation completely disappears, for Fr > 2 (bottom), a roll-wave appears, whose amplitude depends on the initial state (see fig. 8).

Conclusion

This preliminary study of overland flow due to rainfall events clearly enlights several specific difficulties. First, from the numerical point of view, it seems that the apparent topography method, which was designed in order to catch steady states, is not adapted for wet/dry transitions. The semi-implicit treatment seems to be better in the problems we consider and gives good results compared to experimental data. Next, the model itself has to be improved, in particular regarding the empirical friction laws we used, which were not developed in this hydrological context. Finally, more realistic situations require infiltration and two-dimensional simulations, which are in progress and already validated on analytical solutions. This will be again compared with experimental data, as for the flume test.

Figure 1 :

 1 Figure 1: Steady state solution, subcritical inflow and outflow: apparent topography +, semiimplicit ×, analytical -.

Figure 2 :

 2 Figure 2: Left: water front velocities at t = 200 s: apparent topography (+), semi-implicit treatment (×). Right: water front height at t = 200s., semi-implicit

Figure 3 :

 3 Figure 3: Steady state solution, numerical (symbols) vs analytical (lines). Left: subcritical inflow and supercritical outflow, right: supercritical inflow and subcritical outflow.

Figure 4 :

 4 Figure 4: Experimental configuration.

Figure 5 :

 5 Figure 5: Comparison between experimental measures (+) and numerical results (-).

Figure 6 :

 6 Figure 6: Computed rainfall hydrographs for Darcy-Weisbach's law (DW). Left: apparent topography method (AT). Right: semi-implicit scheme.

Figure 7 :

 7 Figure 7: Perturbed initial and final (t = 200s.) states for Fr=1.5 (top left), Fr=2 (top right), Fr=3.7 (bottom).

Figure 8 :

 8 Figure 8: Convergence rate to the final state, Fr≤ 2 (left), Fr>2 (right).

Acknowledgements

This work was partially supported by the ANR grant "METHODE" #ANR-07-BLAN-0232.