Loris Marchal
email: loris.marchal@ens-lyon.fr

Samuel Mccauley

Bertrand Simon
email: bertrand.simon@ens-lyon.fr

Frédéric Vivien
email: frederic.vivien@ens-lyon.fr

Minimizing I/Os in Out-of-Core Task Tree Scheduling

Scientific applications are usually described as directed acyclic graphs, where nodes represent tasks and edges represent dependencies between tasks. For some applications, such as the multifrontal method of sparse matrix factorization, this graph is a tree: each task produces a single output data, used by a single task (its parent in the tree).

We focus on the case when the data manipulated by tasks have a large size, which is especially the case in the multifrontal method. To process a task, both its inputs and its output must fit in the main memory. Moreover, output results of tasks have to be stored between their production and their use by the parent task. It may therefore happen, during an execution, that not all data fit together in memory. In particular, this is the case if the total available memory is smaller than the minimum memory required to process the whole tree. In such a case, some data have to be temporarily written to disk and read afterwards. These Input/Output (I/O) operations are very expensive; hence, the need to minimize them.

We revisit this open problem in this paper. Specifically, our goal is to minimize the total volume of I/O while processing a given task tree. We first formalize and generalize known results, then prove that existing solutions can be arbitrarily worse than optimal. Finally, we propose a novel heuristic algorithm, based on the optimal tree traversal for memory minimization. We demonstrate good performance of this new heuristic through simulations on both synthetic trees and realistic trees built from actual sparse matrices.

I. INTRODUCTION

Parallel workloads are often modeled as task graphs, where nodes represent tasks and edges represent the dependencies between tasks. There is an abundant literature on task graph scheduling when the objective is to minimize the total completion time, or makespan. However, with the increase of the size of the data to be processed, the memory footprint of the application can have a dramatic impact on the algorithm execution time, and thus needs to be optimized. When handling very large data, the available main memory may be too small to simultaneously handle all data needed by the computation. In this case, we have to resort to using disk as a secondary storage, which is sometimes known as out-of-core execution. The cost of the I/O operation to transfer data from and to the disk is known to be several orders of magnitude larger than the cost of accessing the main memory. Thus, in the case of out-of-core execution, it is a natural objective to minimize the total volume of I/O.

In the present paper, we consider the parallel scheduling of rooted in-trees. The vertices of the trees represent computational tasks, and the edges of the trees represent the dependencies between these tasks. Tasks are defined by their input and output data. Each task uses all the data produced by its children to output new data for its parent. In particular, a task must have enough available memory to fit the input from all its children.

The motivation for this work comes from numerical linear algebra, and especially the factorization of sparse matrices using direct multifrontal methods [START_REF] Davis | Direct Methods for Sparse Linear Systems, ser. Fundamentals of Algorithms[END_REF]. During the factorization, the computations are organized as a tree workflow called an elimination tree, and the huge size of the data involved makes it absolutely necessary to reduce the memory requirement of the factorization. Note that we consider here that no numerical pivoting is performed during the factorization, and thus that the structure of the tree, as well as the size of the data, are known before the computation really happens.

It is known that the problem of minimizing the peak memory M peak of a tree traversal, that is, the minimum amount of memory needed to process a tree, is polynomial [START_REF] Liu | An application of generalized tree pebbling to sparse matrix factorization[END_REF], [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF]. However, it may well happen that the available amount of memory M is smaller than the peak memory M peak . In this case, we have to decide which data, or part of data, have to be written to disk. In a previous study [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF], we have focused on the case when the data cannot be partially written to disk, and we proved that this variant of the problem was NP-complete. However, it is usually possible to split data that reside in memory, and write only part of it to the disk if needed. This is for instance what is done using paging: all data are divided in same-size pages, which can be moved from main memory to secondary storage when needed. Since all modern computer systems implement paging, it is natural to consider it when minimizing the I/O volume.

Note that as in [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF], the present study does not directly focus on parallel algorithms. However, parallel processing is the ultimate motivation for this work: complex scientific applications using large data such as multifrontal sparse matrix factorization always make use of parallel platforms. Most involved scheduling schemes combine data parallelism (a task uses multiple processors) and tree parallelism (several tasks are processed in parallel). We indeed have studied such a problem for peak memory minimization [START_REF] Eyraud-Dubois | Parallel scheduling of task trees with limited memory[END_REF]. However, one cannot hope to achieve good results for the minimization of I/O volume in a parallel settings until the sequential problem is well understood, which is not yet the case. The present paper is therefore a step towards understanding the sequential version of this problem.

The main contributions of this work are:

• A formalization in a common framework of the results scattered in the literature; • A proof of the dominance of post-order traversals when trees are homogeneous (all output data have the same size), knowing that an algorithm to compute the best post-order traversal has been proposed by E. Agullo [START_REF] Agullo | On the out-of-core factorization of large sparse matrices[END_REF]. • A proof that neither the best post-order traversal nor the memory-peak minimization algorithms are approximation algorithms for minimizing the I/O volume; • A new heuristic that achieves good performance both on synthetic and actual trees as shown through simulations.

The rest of this paper is organized as follows. We give an overview of the related work in Section II. Then in Section III we formalize our model and present elementary results. Existing solutions are studied in Section IV before a new one is introduced in Section V and evaluated through simulations in Section VI. We finally conclude and present future directions in Section VII.

II. RELATED WORK

Memory and storage have always been a limited parameter for large computations, as outlined by the pioneering work of Sethi and Ullman [START_REF] Sethi | The generation of optimal code for arithmetic expressions[END_REF] on register allocation for task trees. In the realm of sparse direct solvers, the problem of scheduling a tree so as to minimize peak memory has first been investigated by Liu [START_REF] Liu | On the storage requirement in the out-of-core multifrontal method for sparse factorization[END_REF] in the sequential case: he proposed an algorithm to find a peak-memory minimizing traversal of a task tree when the traversal is required to correspond to a postorder traversal of the tree. A postorder traversal requires that each subtree of a given node must be fully processed before the processing of another subtree can begin. A follow-up study [START_REF] Liu | An application of generalized tree pebbling to sparse matrix factorization[END_REF] presents an optimal algorithm to solve the general problem, without the postorder constraint on the traversal. Postorder traversals are known to be arbitrarily worse than optimal traversals for memory minimization [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF]. However, they are very natural and straightforward solutions to this problem, as they allow to fully process one subtree before starting a new one. Therefore, they are widely used in sparse matrix software like MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF], [START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF], and achieve good performance on actual elimination trees [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF]. Note that peak memory minimization is still a crucial question for direct solvers, as highlighted by Agullo et al. [START_REF] Agullo | Robust memory-aware mappings for parallel multifrontal factorizations[END_REF], who study the effect of processor mapping on memory consumption for multifrontal methods.

As mentioned in the introduction, the problem of minimizing the I/O volume when traversing a tree has been studied in [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF] with the constraint that each data either stays in the memory or has to be written wholly to disk. We study here the case when we have the option to store part of the data, which is also the topic of E. Agullo's PhD. thesis [START_REF] Agullo | On the out-of-core factorization of large sparse matrices[END_REF]. In his thesis, Agullo exhibits the best postorder traversal for the minimization of the I/O volume, which we adapt to our model in Section IV-A. He also studies numerous variants of the model that are important for direct solvers, as well as other memory management issues-both for sequential and parallel processing. Based on these preliminaries, he finally presents an out-of-core version of the MUMPS solver.

Finally, out-of-core execution is a well-know approach for computing on large data, especially (but not only) in linear algebra [START_REF] Toledo | A survey of out-of-core algorithms in numerical linear algebra[END_REF], [START_REF] Saule | An out-of-core task-based middleware for data-intensive scientific computing[END_REF].

III. PROBLEM MODELING AND BASIC RESULTS

A. Model and notation

As introduced above, we assume that we have an available memory (or primary storage) of limited size M , and a disk (or secondary storage) of unlimited size.

We consider a workflow of tasks whose precedence constraints are modeled by a tree of tasks G = (V, E ). Its nodes v ∈ V represent tasks and its edges e ∈ E represent dependencies. All dependencies are directed toward the root (denoted by root): a node can only be executed after the termination of all its children. The output data of a node i occupies a size w i in the main memory. This data may be written totally or partially to the disk after task i produces it. In order for a node to be executed, the output data of all its children must be entirely stored in the main memory. An amount of memory m can be moved between the memory and the disk at a cost of m I/O operations, regardless of which data it corresponds to. We assume that all memory values (M , w i ) are given in an appropriate unit (such as kilobytes) and are integers. We divide the main memory into slots, where each slot holds one such unit of memory.

At the beginning of the computation of a task i , the output data of i 's children must be in memory, while at the end of its computation, its own output data must be in memory. The amount of memory needed in order to execute node i is thus wi = max w i , ( j ,i )∈E w j .

We assume that M is at least as large as every wi , as otherwise the tree cannot be processed.

Our objective is to find a solution minimizing the total I/O volume. A solution needs to give the order in which nodes should be executed, and how much of each node should be written out during I/O operations. In particular, for a tree of n tasks, we define a solution to our problem as a permutation σ of [1 . . . n] and a function τ. We call such a solution a traversal. The permutation σ represents the schedule of the nodes, that is, σ(i ) = t means that task i is computed at step t , while the function τ represents the amount of I/O for each data: τ(i ) = m means that m among w i units of the output data of task i are written to disk (then we assume they are written as soon as task i completes). Note that we do not need to clarify which part of the data is written to disk, as our cost function only depends on the volume. Besides, we assume that when τ(i ) = 0, the write operation on the output data of task i is performed right after task i completes (and produces the data), and the read operation is performed just before the use of this data by task i 's parent, as any other I/O scheme would use more memory at some time step for the same I/O volume. Finally, since there are as many read than write operations, we only count the write operations.

In order for a traversal to be valid, it must respect the following conditions:

• Tasks are processed in a topological order:

∀(i , j ) ∈ E , σ(i ) < σ( j );
We say that a node i of parent j is considered active at step t under the schedule σ if σ(i ) < t < σ( j ). This means that its output data is either partially in memory and/or partially written to disk at time t . • The amount of data written to disk never exceeds the size of the data:

∀i ∈ V, 0 ≤ τ(i ) ≤ w i ;
• Enough memory remains available for the processing of each task (taking into account active nodes):

∀i ∈ G, (k,p)∈E σ(k)<σ(i )<σ(p) (w k -τ(k)) ≤ M -wi .
The problem we are considering in this paper, called MINIO, is to find a valid traversal that minimizes the total amount of I/O, given by i ∈G τ(i ).

We formally define a postorder traversal as a traversal σ such that, for any node i and for any node k outside the subtree T i rooted at i , we have either

∀ j ∈ T i , σ(k) < σ( j ) or ∀ j ∈ T i , σ( j ) < σ(k).

B. Towards a compact solution

Although a traversal is described by both the schedule σ and the I/O function τ, the following results show that one can be deduced from the other. The first result is adapted from [5, Property 2.1], which has the same result limited to postorder traversals (see Section II). It states that given a schedule σ, it is easy to derive a I/O scheme τ which minimizes the I/O volume of the traversal (σ, τ). Theorem 1. We consider a tree G, a memory bound M , and a schedule σ. The I/O function τ following the Furthest in the Future policy achieves the best performance under σ.

The I/O function τ following the Furthest in the Future (FiF) policy is defined as follows: during the execution of σ, whenever the memory exceeds the limit M , I/O operations are performed on the active nodes which will remain active the furthest in the future, i.e., whose execution come last in the schedule σ. This result is similar to Belady's rule which states the optimality of the offline MIN cache replacement [START_REF] Belady | A study of replacement algorithms for a virtual-storage computer[END_REF], [START_REF] Lee | A simple proof of optimality for the MIN cache replacement policy[END_REF], that evicts from the cache the data which is used the latest.

Proof. Given a tree G, a memory bound M , a schedule σ and a I/O function τ that does not respect the FiF policy, it is straightforward to transform τ into another I/O function τ following the rule. Consider the first step when an I/O is performed on a data i that is not the last one to be used among active data. Let j denote the last data used among active ones. We can safely increase τ ( j ) and decrease τ (i ) until either τ ( j ) = w j or τ (i ) = 0. As j is active longer than i is, the memory freed by τ is available for a longer time than the one freed by τ, which keeps the traversal valid. Furthermore, by repeating this transformation, we produce an I/O function which respects the FiF policy.

On the other hand, if we have an I/O function τ describing how much of each node is written to disk, we can compute a schedule σ such that (σ, τ) is a valid traversal (if such a schedule exists).

Theorem 2. We consider a tree G, a memory bound M , and an I/O function τ for which there exists a valid schedule. Such a schedule can be computed in polynomial time.

The proof of this result is delegated to Section V where we use a similar method to derive a heuristic: once we know where the I/O operations take place, we may transform the tree by expanding some nodes to make these I/O operations explicit within the tree structure. If a valid traversal using τ exists, the resulting tree may be completely scheduled without any additional I/O, and such a schedule can be computed using an optimal scheduling algorithm for memory minimization.

Both previous results allow us to describe solutions in a more compact format (as either a schedule or an I/O function). However, this does make the problem less combinatorial: there are n! possible schedules and already 2 n τ functions if we restrict only to functions such that τ(i ) = 0 or w i .

C. Related algorithms

As mentioned in Section II, the problem of minimizing the peak memory, denoted MINMEM, is strongly related to our problem, and has been extensively studied. In this problem, the available memory is unbounded (which means no I/Os are required) and we look for a schedule that minimizes the peak memory, i.e., the maximal amount of memory used at any time during the execution. There are at least two important algorithms for this problem, which we use in the present paper:

• It is possible to compute a schedule minimizing the peak-memory in polynomial time, as proved by Liu [START_REF] Liu | An application of generalized tree pebbling to sparse matrix factorization[END_REF].

We denote such an algorithm by OPTMINMEM. • The best postorder traversal for peak-memory minimization can also be computed in polynomial time [START_REF] Liu | On the storage requirement in the out-of-core multifrontal method for sparse factorization[END_REF].

We will refer to this algorithm by POSTORDERMINMEM.

IV. EXISTING SOLUTIONS ARE NOT SATISFACTORY

We now detail two existing solutions for the MINIO problem. The first one is the best postorder traversal proposed by Agullo [START_REF] Agullo | On the out-of-core factorization of large sparse matrices[END_REF]. The second consists in using the optimal traversal for MINMEM proposed by Liu [START_REF] Liu | An application of generalized tree pebbling to sparse matrix factorization[END_REF] and then to apply Theorem 1 to obtain a valid traversal. After presenting these algorithms, we prove that none of them has a constant competitive factor compared to the optimal traversal.

A. Computing the best postorder traversal

For the sake of completeness, we present the adaption to our model of the algorithm computing the best postorder traversal for MINIO from [START_REF] Agullo | On the out-of-core factorization of large sparse matrices[END_REF]. Recall that in a postorder traversal, when a node is processed, its whole subtree must be processed before any other external node may be started. Given a node i and a postorder schedule σ, we first recursively define S i as the storage requirement of the subtree T i rooted at i . Let Chil(i ) be the children of i . Then:

S i = max     w i , max j ∈Chil(i )     S j + k∈Chil(i ) σ(k)<σ( j ) w k         .
This expression represents the maximum memory peak reached during the execution. If the peak is obtained at the end of the execution, it is then equal to w i . Otherwise, it appears during the execution of the subtree of some child j . In this case, the peak is composed of the weights of the children already processed, plus the peak S j of T j .

We may now consider A i = min(M , S i ), which represents the amount of main memory used for the out-of-core execution of the subtree T i by σ. We recursively define V i as the volume of I/Os performed by σ during the execution T i when I/O operations are done using the FiF policy:

V i = max 0, max j ∈Chil(i ) A j + σ(k)<σ( j ) k∈Chil(i ) w k -M + j ∈Chil(i ) V j .
The expression of V i has a similar structure to the expression of S i . No I/Os can be incurred when only the root i is in memory, hence w i has no effect here. The second term accounts for the I/Os incurred on the children of i . Indeed, during the execution of node j , some parts of children of i must be written to disk if the memory peak exceeds M , and this quantity is at least A j + σ(k)<σ( j ) k∈Chil(i ) w k -M . The last term accounts for the I/Os occurring inside the subtrees. Note that such I/Os can only happen if the memory peak of the subtree exceeds M .

It remains to determine which postorder traversal minimizes the quantity V root . Note that the only term sensitive to the ordering of the children of i in the expression of V i is max

j ∈Chil(i ) A j + σ(k)<σ( j ) k∈Chil(i )
w k . Theorem 3 states that sorting the children of i by decreasing values of A jw j achieves the minimum V i . Theorem 3 (Lemma 3.1 in [START_REF] Liu | On the storage requirement in the out-of-core multifrontal method for sparse factorization[END_REF]). Given a set of values (x i , y i ) 1≤i ≤n , the minimum value of max 1≤i ≤n x i + i -1 j =1 y j is obtained by sorting the sequence (x i , y i ) in decreasing order of x iy i .

Therefore, the postorder traversal that processes the children nodes by decreasing order of A iw i minimizes the I/O cost among all postorder traversals. This traversal is described in Algorithm 1, initially called with r = root, and will be referred to as POSTORDERMINIO. Note that in the algorithm ⊕ refers to the concatenation operation on lists.

Algorithm 1: POSTORDERMINIO (G, r )

Output: a tree G and a node r in G Output: an ordered list r of the nodes in the subtree rooted at r , corresponding to a postorder

1 foreach i child of r do 2 i ← POSTORDERMINIO(G, i ) 3 Compute the A i value using postorder i 4 r ← 5 for i child of r by decreasing value of A i -w i do 6 r ← r ⊕ i 7 r ← r ⊕ {r } 8 return r

B. POSTORDERMINIO is optimal on homogeneous trees

In this subsection we focus on homogeneous trees, which are defined as trees in which all nodes have output data of size one. We will show that POSTORDERMINIO is optimal on these homogeneous trees, i.e., that it performs the minimum number of I/Os. This generalizes a result of Sethi and Ullman [START_REF] Sethi | The generation of optimal code for arithmetic expressions[END_REF] for homogeneous binary trees. Theorem 4. POSTORDERMINIO is optimal for homogeneous trees.

In order to prove this theorem, we need first to define some labels on the nodes of a tree. Let T be any homogeneous tree (w v = 1 for all nodes v of T ). In the following definitions, whenever v is a node of T with k children, v 1 , . . . , v k will be its children.

Memory bound l (v).

For each node v of T , we recursively define a label l (v) which represents the minimum amount of memory necessary to execute the subtree T (v) rooted at v without performing any I/Os:

l (v) =        0 if v is a leaf max 1≤i ≤k (l (v i ) + i -1)
otherwise and ordering the children such that l

(v i ) ≥ l (v i +1 ) for 1 ≤ i ≤ k -1
We call POSTORDER one postorder schedule that executes the children of any node by non-increasing llabels (ties being arbitrarily broken). Intuitively, under POSTORDER, while computing the i -th child, we have i -1 extra nodes in memory, each of size one, so we need l (v i ) + (i -1) memory slots in total.

I/O indicator c(v). If v i is a child of v, intuitively, c(v i )
represents the number of children of v written to disk by POSTORDER during the execution of T (v i ). This number can be either 0 or 1. We set c(v 1 ) = 0 and

c(v i ) = 0 if l (v i ) + 1≤ j <i (1 -c(v j )) ≤ M 1 otherwise.
We set c(root) = 0. To ease the writing of some proofs, we use the notation

m(v i ) = 1≤ j <i (1 -c(v j )).
Thus m(v i ) represents the number of children of v in memory right before v i is executed. Note that m(v 1 ) = 0 and m

(v i ) = (1-c(v 1 ))+ 2≤ j <i (1-c(v j )) ≥ (1-c(v 1 )) = 1 for 2 ≤ i ≤ k. I/O volumes w(v) and W (T (v)). w(v) represents the total
number of children of v stored by POSTORDER:

w(v) = k i =1 c(v i ) = k i =2 c(v i ).
Finally, for a given node v, we define W (T (v)) on the subtree rooted at v:

W (T (v)) = c(v) + µ∈T (v)
w(µ).

W (T (v)) intuitively represents the total volume of communications performed during the execution of the tree T (v) by POSTORDER.

We first state the correctness of the l -labels and the optimality of POSTORDER for the MINMEM problem. Lemma 1. With infinite memory, POSTORDER uses l (n) slots to compute the subtree rooted at node n.

Proof. This result follows from the definition of the labels l (v).

Lemma 2. With infinite memory, any schedule uses at least l (v) slots to compute the subtree rooted at v.

Proof. We prove this result by induction on the size of T (v). If v is a leaf, the result holds (l (v) = 1).

Otherwise, we assume the lemma to be true for the subtrees rooted at the children v 1 , . . . , v k of v. We consider the schedule returned by MINMEM. The memory peak inherent to the execution of a subtree T (v i ) is equal to l (v i ) by the induction hypothesis. Assume without loss of generality that the children of v are ordered such that MINMEM first computes a node of T (v 1 ), then the next executed node not in T (v 1 ) is in T (v 2 ), then the next executed node neither in T (v 1 ) nor in T (v 2 ) is in T (v 3 ), and so on. Then, the memory peak reached during the execution of T (v i ) is at least l (v i ) + (i -1) because, in addition to T (v i ), at least i -1 subtrees have been partially executed: T (v 1 ), ..., T (v i -1 ). Finally, the total memory peak is at least equal to max 1≤i ≤k (l (v i )+i -1). By Theorem 3, this quantity is minimized when the nodes are ordered by nonincreasing values of l (v i ). Hence, the total memory peak is at least l (v).

We now state the performance of POSTORDER for the MINIO problem (I/Os are done following the FiF policy).

Lemma 3.

With bounded memory M , POSTORDER computes a given tree T using W (T ) I/Os.

Proof. We prove this result by induction on the size of T . We introduce new notation: for any node v of T we define

W (v) as W (v) = W (T (v)) -c(v).
In other words, W (v) intuitively represents the total volume of communications performed during the execution of the tree T (v) if we had nothing to execute but T (v) (in practice T (v) may be a strict subtree of T and, therefore, the execution of T (v) in the midst of the execution of T can induce more communications). Note that W (v) = W (T (v)) if v is the root of T . We prove by induction on the size of T (v) that W (v) I/Os are performed during the execution of T (v).

Let us assume that v is a leaf. Then W (v) = 0. Because we have assumed (in Section III-A) that M was large enough for a single node to be processed without I/Os, c(v) = 0 and thus

W (T (v)) = 0 = W (v) + c(v).
Now assume that v is not a leaf. By the induction hypothesis, for any i ∈ [1; k], POSTORDER executes the tree T (v i ) alone using W (v i ) I/Os. We prove that to process the tree T (v i ), after the trees T (v 1 ) through T (v i -1 ) were processed, we need to perform

W (T (v i )) = W (v i ) + c(v i ) I/Os. Let us consider the (i + 1)-th child of v. If c(v i +1 ) = 0, then l (v i +1 ) + 1≤ j <i +1 (1 -c(v j )) ≤ M .
Then, according to Lemma 1, no I/Os are required to execute T (v i +1 ) under POSTORDER even after the processing of T (v 1 ) through T (v i ). Indeed, before the start of the processing of T (v i +1 ) the memory contains exactly 1≤

j <i +1 (1 -c(v j )) nodes. Therefore W (v i +1 ) = c(v i +1 ) = W (T (v i +1 )) = 0.
We are now in the case c(

v i +1 ) = 1; thus l (v i +1 ) + 1≤ j <i +1 (1 -c(v j )) > M . Recall that for l ∈ [1; i ], l (v l ) ≥ l (v i +1 ). Thus, if l (v i +1 ) ≥ M , then for l ∈ [2; i ], l (v l ) ≥ M and c(v l ) = 1 (because m(v l ) ≥ (1 -c(v 1 )) = 1)
. Therefore, after the completion of T (v i ) there is only one node remaining in the memory: v i . Then with a single I/O POSTORDER writes v i to disk, the memory is empty and T (v i +1 ) can then be processed with W (v i +1 ) I/Os, giving a total of W (v i +1 ) + c(v i +1 ) = W (T (v i +1 )) I/Os. The only remaining case is the case l (v i +1 ) < M . The processing of T (v i ) requires at least l (v i +1 ) empty memory slots because l (v i ) ≥ l (v i +1 ). Hence, after the completion of T (v i ) there are at least l (v i +1 ) -1 empty memory slots (the memory including the node v i itself ). Then with a single I/O POSTORDER writes v i to disk and there are enough empty memory slots to process T (v i +1 ) without any additional I/Os. Therefore W (T (v i +1 )) = 1 = W (v i +1 ) + c(v i +1 ). This concludes the proof. Lemma 5 relies on the following intermediate result.

Lemma 4. Consider a node v of a tree T with a child, a, whose label l (a) satisfies l (a) > M . Now, consider any tree T identical to T , except that the subtree rooted at a has been replaced by any tree whose new label l (a) satisfies l (a) ≤ l (a) and l (a) ≥ M . Then w (v) = w(v).

Proof. Let v 1 , . . . , v k be the children of v, ordered so that l (v 1 ) ≥ • • • ≥ l (v k ). Let j be the index of a: a = v j . As the label of a in T , l (a), is not larger than l (a), we can have l (a) < l (v j +1 ). Therefore, we define another ordering of the children of v denoted by v 1 , . . . , v k such that l (v 1 ) ≥ • • • ≥ l (v k ). Let j be the index of a in this ordering:

v j = a = v j .
Note that j ≥ j . For i ∈ [ j + 1; j ], we have v i = v i -1 ; at j , we have v j = v j ; and for i ∉ [ j ; j ], we have v i = v i .

If j = 1 then j = 1. This case means that a remains the node with the largest label. The labels of the other children of v remain unchanged. Because c(v 1 ) = c (v 1 ) = 0 by definition, then c (v i ) = c(v i ) for any child v i of v and, thus, w(v) is equal to w (v).

Let us now consider the case j > 1. From what precedes,

v j -1 = v j . Then l (v j ) = l (v j -1 ) ≥ l (v j ) = l (a) ≥ M . However, for any i ∈ [1; j ], l (v i ) ≥ l (v j ) ≥ M and l (v i ) ≥ l (v j ) ≥ M . Therefore, for any i ∈ [2; j ], l (v i ) + m (v i ) > M (because m (v i ) ≥ 1 -c (v 1 ) = 1 ) and, thus, c (v i ) = 1. Similarly, for any i ∈ [2; j ], l (v i ) + m(v i ) > M (because m(v i ) ≥ m(v 1 ) = 1 ) and, thus, c(v i ) = 1. Therefore, for i ∈ [1; j ], c(v i ) = c (v i ). Then, for i ∈ [ j + 1; k], v i = v i , m(v i ) = m (v i ), and c(v i ) = c (v i ) by an obvious induction. Therefore, w (v) = k i =2 c (v i ) = k i =2 c(v i ) = w(v).
The following lemma gives a lower bound on the I/Os performed by any schedule.

Lemma 5. No schedule can compute a tree T performing strictly less than W (T ) I/Os.

Proof. This proof has been deferred to the companion research report [START_REF] Marchal | Minimizing i/os in out-of-core task tree scheduling[END_REF], due to its technicality and its length. The main ideas of this proof are thus summarized here. The result is proved by induction on the size of the tree. The case where no I/O is required is deduced from Lemma 2. Now, we consider a tree T for which any schedule performs at least one I/O, and an optimal schedule P on this tree. We focus on the first node s to be stored under this schedule, and define the tree T in which T (s) is replaced by s. Using the induction hypothesis, we know that any schedule on T , including the restriction of P on T , performs at least W (T ) I/Os. Therefore, we deduce that P performs at least W (T )+1 I/Os on T . Thus, the objective is to prove that W (T ) ≥ W (T ) -1.

In order to prove this result, we focus on the closest ancestor of s to have a label l larger than M , and denote it as µ. We first prove that in the new tree T , we have l (µ) ≥ M . This means, by Lemma 4, that the w labels of the ancestors of µ are unchanged in T . Then, we prove through an extensive case study that w(µ) in T cannot be smaller than w(µ) in T minus one. Finally, we conclude that all the others w labels are equal in T and in T , therefore, W (T ) ≥ W (T ) -1.

We are now ready to prove Theorem 4.

Proof of Theorem 4. Because of Lemma 3 and of Lemma 5, POSTORDER is optimal for homogeneous trees. However, POSTORDERMINIO is a post-order that minimizes the volume of I/O operations. Hence, it is also optimal for homogeneous trees.

Note that, on homogeneous trees, POSTORDER and POS-TORDERMINIO are almost identical: POSTORDER sorts children by non-increasing l i , while POSTORDERMINIO sorts them by non-increasing A i = mi n(M , l i -1). In particular, for children with l i > M , the order is not significant for POSTORDERMINIO.

C. Postorder traversals are not competitive

Previous research has shown that the best postorder traversal for the MINMEM problem is arbitrarily far from the optimal traversal [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF]. We prove here that postorder traversals may also have bad performance for the MINIO problem. More specifically, we prove that there exist problem instances on which POSTORDERMINIO performs arbitrarily more I/O than the optimal I/O amount. We could exhibit an example where the optimal traversal does not perform any I/O and POSTORDERMINIO performs some I/O, but we rather present a more general example where the optimal traversal performs some I/O: in the following example, the optimal traversal requires 1 I/O, when POSTORDERMINIO requires Ω(nM ) I/Os. The tree used in this instance is depicted on Figure 1(a).

It is possible to traverse the tree of Figure 1(a) with a memory of size M using only a single I/O, by executing the nodes in increasing order of the labels next to the nodes. After processing the minimal subtree including the two leftmost leaves, our strategy is to process leaves from left to right. Before processing a new leaf, we complete the previous subtree up to a node of weight 1; this way the leaf and the actives nodes can both fit in memory.

On the other hand, the best postorder traversal must perform a volume of I/O equal to M /2 -1 before it can start any leaf except for the first leaf it processes. This is because the least common ancestor of any two leaf nodes has two children of size M /2, and all leaves have size at least M -1. Thus, any postorder traversal incurs at least M /2 -1 I/Os for all but one leaf node (3M /2 -2 for the example here). We can extend the tree in Figure 1(a): we replace root by a node of size 1, add to it a parent of size

root 15 M /2 14 1 11 M /2 10 1 7 M /2 6 1 2 M 1 M /2 5 1 4 M 3 M /2 9 M -1 8 M /2 13 M -1 12
(a) Example of a tree showing that POSTORDERMINIO is not an approximation algorithm.

root 2k 4k + 4 3k 4k + 3 2k -1 4k -2 3k + 1 4k -3 • • • k 4 4k 3 2k 4k + 2 3k 4k + 1 2k -1 4k 3k + 1 4k -1 • • • k 2 4k 1
(c) Example of a tree showing that OPTMINMEM is not an approximation algorithm (M = 4k). M /2 which is the left child of the new root; the right child of the new root is then a chain containing a leaf of size M -1 and its parent of size M /2. Doing this repeatedly until n nodes are used gives the lower bound of Ω(nM ). Therefore, POSTORDERMINIO is not constant-factor competitive.

D. OPTMINMEM is not competitive

Minimizing the amount of I/O in an out-of-core execution seems close to minimizing the peak memory when the memory in unbounded. Thus, in order to derive a good solution for MINIO, it seems reasonable to use an optimal algorithm for MINMEM, such as the OPTMINMEM algorithm presented by Liu [START_REF] Liu | An application of generalized tree pebbling to sparse matrix factorization[END_REF], to compute a schedule σ and then to perform I/Os using the FiF policy. In the following, we also use OPTMINMEM to denote this strategy for MINIO. We prove here that there exist problem instances on which this strategy will also perform arbitrarily more I/Os than the optimal traversal. We first exhibit in Figure 1(b) a tree showing that OPT-MINMEM does not always lead to minimum I/Os in our model. Let M = 6. The tree of Figure 1(b) can be completed with 3 I/Os, by doing one chain after the other. This corresponds to a peak memory of 9. But OPTMINMEM achieves a peak memory of 8 at the cost of 4 I/Os by executing the nodes in increasing order of the labels next to the nodes.

This example can be extended to show that OPTMIN-MEM may perform arbitrarily more I/Os than the optimal strategy. The extended tree is illustrated on Figure 1(c). It contains two identical chains of length 2k + 2, for a given parameter k, and the memory size is set to 4k. The weights of the tasks in each chain (in order from root to leaf) are defined by interleaving two sequences: {2k, 2k -1, . . . , k} and {3k, 3k + 1, . . . , 4k}. As above, it is possible to schedule this tree with only 2k I/Os, but with a memory peak of 6k, by computing first one entire chain. However, OPTMINMEM achieves a memory peak of 5k by switching chains on each node with a weight smaller than 2k, as represented by the labels besides the nodes. Doing so, OPTMINMEM incurs k I/Os on each of the k +1 smallest nodes, leading to a cost of k(k +1) I/Os. The competitive ratio is then larger than k/2. Therefore, OPTMINMEM is not constant-factor competitive in the MINIO problem.

E. Complexity unknown

As shown above, polynomial-time approaches based on similar problems fail to even give a constant-competitive ratio. The main issue facing a polynomial approach is the highly nonlocal aspect of the optimal solution. For example, since postorder traversals are not optimal, it may be highly useful to stop at intermediate points of a subtree's execution in order to process entirely different subtrees.

We conjecture that this problem is NP-hard due to these difficult dependencies. As mentioned above, if we require entire nodes to be written to disk, the problem has been shown to be NP-hard by reduction to Partition [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF]. However, this proof highly depends on indivisible nodes, rather than on the tree's recursive structure. Taking advantage of the structure of our problem to give an NP-hardness result could lead to an interesting understanding of optimal solutions, and possibly further heuristics. We leave this as an open problem.

w i ⇒ w i w i -τ(i ) w i

V. HEURISTIC

We now move to the design of a novel heuristic FULL-RECEXPAND whose goal is to improve the performance of OPTMINMEM for the MINIO problem. The main idea of this heuristic is to run OPTMINMEM several times: when we detect that some I/O is needed on some node, we force this I/O by transforming the tree. This way, the following iterations of OPTMINMEM will benefit from the knowledge of this I/O. We continue transforming the tree until no more I/Os are necessary.

In order to enforce I/Os, we use the technique of expanding a node (illustrated on Figure 2). Under an I/O function τ, we define the expansion of a node i as the substitution of this node by a chain of three nodes i 1 , i 2 , i 3 of respective weights w i , w i -τ(i ) and w i . The expansion of a node actually mimics the action of executing I/Os: the weight of the three tasks represent which amount of main memory is occupied by this node 1) when it is first completed (w i 1 = w i ), 2) when part of it is moved to disk (w i 2 = w i -τ(i )), and 3) when the whole data is transferred back to main memory (w i 3 = w i ).

This technique first allows us to prove Theorem 2, which states that given an I/O function τ, we can find a schedule σ such that (σ, τ) is a valid traversal if there exists one.

Proof of Theorem 2. Consider the tree G obtained from G by expanding all the nodes for which τ is not null. Then, consider the schedule σ obtained by OPTMINMEM on G , and let σ be the corresponding schedule on G. Then, the memory used by σ on G during the execution of a node i is the same as the one used by σ on G during the execution of the same node i , or of i 1 if i is expanded. Then, as OPTMINMEM achieves the optimal memory peak on G , we know that σ uses as little main memory as possible under the I/O function τ. Then, (σ, τ) is a valid traversal of G.

The heuristic FULLRECEXPAND is described in Algorithm 2. The main idea of the heuristic is to expand nodes in order to obtain a tree that can be scheduled without I/O, which is equivalent to building an I/O function.

First, the heuristic recursively calls itself on the subtrees rooted at the children of the root, so that each subtree can be scheduled without I/O (but using expansions). Then, the algorithm computes OPTMINMEM on this new tree, and if I/Os are necessary, it determines which node should be expanded next. This selection is the only part where FULLRECEXPAND can deviate from an optimal strategy. Our choice is to select a node on which the FiF policy would incur I/Os; if there are several such nodes, we choose the one whose parent is scheduled the latest. After the expansion, the algorithm recomputes OPTMINMEM on the modified tree, and proceeds until no more I/O are necessary.

At the end of the computation, the returned schedule is obtained by running OPTMINMEM on the final tree computed by FULLRECEXPAND, and by transposing it on the original tree. The I/O performance of this schedule is then equal to the sum of the expansions. FULLRECEXPAND is only a heuristic: it may give suboptimal results but also achieve better performance than OPTMINMEM, as illustrated on several examples in the companion research report [START_REF] Marchal | Minimizing i/os in out-of-core task tree scheduling[END_REF] .

Unfortunately, the complexity of FULLRECEXPAND is not polynomial, as the number of iterations of the while loop at Line 3 cannot be bounded by the number of nodes, but may depend also on their weights. We therefore propose a simpler variant, named RECEXPAND, where the while loop at Line 3 is exited after 2 iterations. In this variant, the resulting tree G might need I/Os to be executed. The final schedule is computed as in FULLRECEXPAND, by running OPTMINMEM on this tree G. We later show that this variant gives results which are very similar to the original version.

VI. NUMERICAL RESULTS

In this section, we compare the performance of the two existing strategies OPTMINMEM and POSTORDERMINIO, and the two proposed heuristics FULLRECEXPAND and RE-CEXPAND. All algorithms are compared through simulations on two datasets described below. Because of its high computational complexity, FULLRECEXPAND is only tested on the first smaller dataset.

A. Datasets

The first dataset, named SYNTH, is composed of 330 instances of synthetic binary trees of 3000 nodes, generated uniformly at random among all binary trees. As we considered small trees, we simply used half-Catalan numbers in order to draw a tree, similarly to the method described at the beginning of [START_REF] Mäkinen | Generating random binary trees-a survey[END_REF]. The memory weight of each task is uniformly drawn from [1; 100]. The second dataset, named TREES, is composed of 329 elimination trees of actual sparse matrices from the University of Florida Sparse Matrix Collection 1 (see [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF] for more details on elimination trees and the data set). Our datasets corresponds to the 329 smallest of the 640 trees presented in [START_REF] Jacquelin | On optimal tree traversals for sparse matrix factorization[END_REF], with trees ranging from 2000 to 40000 nodes.

For each tree of the two datasets, we first computed the minimal memory size necessary to process the tree nodes: LB = max i wi . We also computed the minimal peak memory for an incore execution Peak incore (using OPTMINMEM). We eliminated some trees from the TREES dataset where Peak incore = LB, leaving us with 133 remaining trees in this dataset. In all other cases, note that the possible range for the memory bound M such that some I/Os are necessary is [LB, Peak incore -1]. We chose to set M to the middle of this interval M = (LB + Peak incore -1)/2. In the companion research report [START_REF] Marchal | Minimizing i/os in out-of-core task tree scheduling[END_REF], we relate simulations using other values of M in this interval which show similar trends.

B. Results

Our objective in this study is to minimize the total amount of I/Os needed to process the tree. In order to summarize and compare the performance of the different strategies we choose here to consider the number of I/Os and the memory bound M : performing 10 I/Os when the optimal only needs 1 does not have the same significance if the main memory consists of M = 10 slots or of M = 1000 slots. Therefore, in this section, if a schedule performs k I/Os, we define its performance as (M + k)/M . Then, a schedule with no I/O operations has a performance of 1 while a schedule needing M I/Os has a performance of 2.

In order to compare the performance of these algorithms, we use a generic tool called performance profile [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF]. For a given dataset, we compute the performance of each algorithm on each tree and for each memory limit. Then, instead of computing an average above all the cases, a performance profile reports a cumulative distribution function. Given a heuristic and a threshold (or overhead) τ expressed in percentage, we compute the fraction of test cases in which the performance of this heuristic is at most τ% larger than the best observed performance, and plot these results. Therefore, the higher the curve, the better the method: for instance, for an overhead τ = 5%, the performance profile shows how often a given method lies within 5% of the smallest performance obtained.

The left plot of Figure 3 presents the performance profile of the four heuristics for the complete dataset SYNTH. The first result is the poor performance of POSTORDERMINIO in this dataset: it almost always has at least 50% of overhead, and even a 100% overhead in 75% of the cases. Then, RE-CEXPAND performs far better than OPTMINMEM. The right plot of the figure presents the performance profiles of exclusively OPTMINMEM, RECEXPAND and FULLRECEXPAND. RECEXPAND performs strictly less I/Os than OPTMINMEM on 90% of the instances, and on half of them, OPTMINMEM has a 4% overhead. We can also note that FULLRECEXPAND performs only slightly better than RECEXPAND, but both heuristics are far ahead of OPTMINMEM, so the gain in the complexity of the algorithm is only balanced by a small loss of performance. For instance, RECEXPAND has an overhead of more than 2% over FULLRECEXPAND on only 3% of the instances.

The left plot of Figure 4 presents the performance profiles of the three heuristics POSTORDERMINIO, RECEXPAND and OPTMINMEM for the complete dataset TREES. The first remark is that the three heuristics are equal on more than 90% of the 329 instances. Therefore, we now focus on the right plot, which presents the same performance profile for the 25 cases where the heuristics do not all give equal performance. We can see that the hierarchy is the same as in the previous dataset (RECEXPAND is never outperformed, and OPTMINMEM performs better than POSTORDERMINIO) but with smaller discrepancies between the heuristics. POS-TORDERMINIO and OPTMINMEM respectively have more than 5% of overhead on only 40% and 10% of these instances.

VII. CONCLUSION

In this paper, we revisited the problem of minimizing I/O operations in the out-of-core execution of task trees. We proved that existing solutions allow us to optimally solve the problem when all output data have identical size, but that none of them has a constant competitive factor compared to the optimal solution. We proposed a novel heuristic solution that improves on an existing strategy and is very efficient in practice. Despite our efforts, the complexity of the problem remains open. Determining this complexity would definitely be a major step, although our findings already lay the basis for more advanced studies. These include moving to parallel out-of-core execution (as was already done for parallel incore execution [START_REF] Eyraud-Dubois | Parallel scheduling of task trees with limited memory[END_REF]) as well as designing competitive algorithms for the sequential problem.

VIII. ACKNOWLEDGMENT

  Example of a tree where OPTMINMEM is not optimal for MINIO (M = 6).

Figure 1 :

 1 Figure 1: The label inside node i represents w i . The label next to the nodes indicate in (a) the optimal schedule, and in (b) and (c) the OPTMINMEM schedule.

Figure 2 :

 2 Figure 2: Example of node expansion.

Algorithm 2 :1 foreach child i of r do 2 G 5 τ 6 i

 2256 FULLRECEXPAND (G, r, M ) Input: tree G, root of exploration r Output: Return a tree G r which can be executed without I/O, obtained from G by expanding several nodes i ← FULLRECEXPAND(G, i , M ) 3 G r ← tree formed by the root r and the subtrees G i 4 while OPTMINMEM(G r , r ) needs more than a memory M do ← I/O function obtained from OPTMINMEM(G r , r ) using the FiF policy ← node for which τ(i ) > 0 whose parent is scheduled the latest in OPTMINMEM(G r , r ) 7 modify G r by expanding node i according to τ(i ) 8 return G r

Figure 4 :

 4 Figure 4: Performance profiles for the complete TREES dataset (left) and restricted to instances where the heuristics differ (right)

  1 http://www.cise.ufl.edu/research/sparse/matrices/

This material is based upon research supported by the SOLHAR project operated by the French National Research Agency (ANR) and the Chateaubriand Fellowship of the Office for Science and Technology of the Embassy of France in the United States.


Alexis Maineult

Laurence Jouniaux

Yves Bernabé

Influence of the mineralogical composition on the self-potential response to advection of KCl concentration fronts through sand

published or not. The documents may come

Influence of the mineralogical composition on the self-potential response to advection of KCl concentration fronts through sand

Alexis Maineult, Laurence Jouniaux, Yves Bernabé

INTRODUCTION

Self-potential (SP) monitoring is a geophysical method based on the measurement of voltage fields naturally occurring in the Earth's subsurface. Owing to the sensitivity of the self-potential to variations in groundwater flow, chemistry or temperature [e.g., [START_REF] Corwin | The self-potential method in geothermal exploration[END_REF], the method has been frequently used recently in subsurface studies [e.g., [START_REF] Sandberg | An integrated geophysical investigation of the hydrogeology of an anisotropic unconfined aquifer[END_REF][START_REF] Jardani | Self-potential tomography applied to the determination of cavities[END_REF]. However, the simultaneous action of many different SP sources can make the interpretation difficult in terms of fluid velocity and composition. To better characterize these processes, sand-box experiments involving reproducible, controlled sources have been conducted in the laboratory [START_REF] Ahmad | A laboratory study of streaming potentials[END_REF][START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF][START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF][START_REF] Maineult | Detection of advected, reacting redox fronts from self-potential measurements[END_REF][START_REF] Suski | A sandbox experiment of self-potential signals associated with a pumping test[END_REF][START_REF] Naudet | A sandbox experiment to investigate bacteria-mediated redox processes on self-potential signals[END_REF]. In previous works focusing on the SP response to advective transport of NaCl and FeCl 2 concentration fronts [START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF][START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF][START_REF] Maineult | Detection of advected, reacting redox fronts from self-potential measurements[END_REF], we found that, to first order, the SP signals could be modeled as the sum of: 1) an electrokinetic term, related to local changes of the electrokinetic coupling coefficient with varying fluid concentration, and 2) a fluid junction potential term, caused by the presence of a concentration gradient. Noting that the second term should vanish when ionic species with identical mobilities are used, we tried to test this prediction by running advective transport experiments using a KCl solution (the mobility of K + is equal to 96 % of the mobility of Cl -). However the observations appeared to contradict this theoretical prediction. This unexpected discrepancy led us to perform a new set of KCl advection experiments using a 100 % pure quartz sand instead of the previously employed sand, which contained small amounts of feldspars, clays and micas. The results reported here suggest that the anomalous SP signals can be attributed to chemical reactions between K + and the minor minerals.

3

BACKGROUND

Fluid flow through a porous medium produces a downstream motion of the counter-ions in excess near the charged pore-walls, resulting in a net charge separation and, consequently, in the generation of the so-called electrokinetic field. Assuming no surface conductivity, the electric and hydraulic potentials U e (in V) and H (in m) are related to each other by the Helmholtz-Smoluchowski equation:

* e g U L H H ρ εζ ησ = = ∇ ∇ ∇ (1)
where, ρ, ε, η and σ are the density (kg m -3 ), dielectric constant (F m -1 ), viscosity (Pa s)

and electrical conductivity (S m -1 ) of the fluid, g the gravity acceleration (m s -2 ) and L * the electrokinetic coupling coefficient (V m -1 ). The so-called "zeta-potential" ζ (V) depends on the geochemical properties of the rock, and on the fluid composition and concentration [e.g., [START_REF] Ogilvy | Geophysical studies of water leakages from reservoirs[END_REF][START_REF] Ishido | Experimental and theoretical basis of electrokinetic phenomena in rock-water systems and its applications to geophysics[END_REF][START_REF] Morgan | Streaming potential properties of Westerly granite with applications[END_REF][START_REF] Lorne | Streaming potential measurements 1. Properties of the electrical double layer from crushed rocks[END_REF][START_REF] Guichet | Effect of secondary minerals on electrokinetic phenomena during water-rock interaction[END_REF][START_REF] Hase | ζ potential measurements of volcanic rocks from Aso caldera[END_REF]. Equation (1) implies that any change in the fluid concentration can modify the electrokinetic potential.

As ionic species diffuse along a concentration gradient, a net charge separation occurs if the anion and cation mobilities are different, thus generating a counteracting electrical field. The junction potential U j (V) obeys the Planck-Henderson equation:

( ) R Ae j T U σ σ ϕ σ σ - + - + - = + ∇ ∇ (2)
where σ -and σ + are the anionic and cationic electrical conductivities (such as σ -+ σ + = σ), ϕ the porosity, T the absolute temperature (K), R the molar gas constant (8.314 J mol - 4 (≈1.602 10 -19 C). In the case of a symmetric, monovalent, binary salt such as NaCl or KCl, equation (2) can be simplified as:

* R Ae j u u C T C U C u u C α ϕ - + - + - = = + ∇ ∇ ∇ (3)
where α * is the junction coupling coefficient in a porous medium (V), C the salt concentration (mol m -3 ), and u -and u + the ionic mobilities (m 2 s -1 V -1 ) of anions and cations, respectively.

MATERIALS AND METHOD

The flow-through sand-box

We performed the experiments reported here in the sand-box described in detail in [START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF][START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF] (Figure 1). Briefly, a rectangular container (44.25 × 23.75 × 26.50 cm) was divided into three hydraulically connected compartments. The endreservoirs (5.7-and 6.1-cm long) were used to control the flow, while the 31-cm long inner one was filled up with a 21-cm height water-saturated sand-body. Precautions were taken to ensure packing homogeneity and avoid air trapping.

A uniform, one-dimensional flow field was produced by tilting the container by an angle of 4.4° (corresponding to a mean hydraulic gradient ∇H = -7.7 %) while maintaining the water level in each end-reservoir constant at 20 cm from the box bottom.

Note that about 60 liters of primary fluid were circulated through the sand in a closed circuit for several days before the experiment, to reach steady-state flow and chemical equilibrium conditions. After establishing steady-state flow (i.e. stabilized flowrate, conductivity and electrical potential differences), we changed the setting to open-circuit and generated a sharp concentration pulse at time t = 0 by quickly adding and stirring 6 cm 3 of saturated KCl solution in the upstream reservoir. This volume was small enough to minimize the perturbation of flow (the resulting transitory, upstream water level increase was only ≈0.5 mm). Generating the pulse took less than 10 s. We then stirred the upstream reservoir every 2 minutes to ensure a homogeneous upstream solution. This procedure produced a sharp concentration front followed by an exponentially decaying tail (see Figure 2).

We measured the electrical potential differences (or EPDs) between a reference electrode placed in the upstream reservoir and four electrodes inserted at half height in the sand body, along a median, longitudinal line, at distances y equal to 5, 12, 19 and 26 cm from the upstream reservoir (Figure 1). We used the small-size (i.e., 4 mm in diameter), copper-copper sulfate, unpolarisable electrodes described in [START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF]. To remove the intrinsic potential of each electrode and the constant electrokinetic signal produced by the steady-state flow of the primary fluid, the EPD signals were reduced by their individual initial value (all curves shown hereafter thus represent variations with respect to the initial steady-state flow of the primary fluid).

Sands

The S4M sand from Haguenau, France, contains 96 % quartz, about 3 % potassic feldspars (microcline and sanidine) and slightly more than 1 % micas and clay minerals (essentially altered biotite, muscovite and illite). For the non-clay minerals, the grain diameter ranges in the interval 200 -400 µm, with a mean of 292 µm. The porosity ϕ of several similarly packed sand samples was determined by dry and water-saturated weighting and ranged between 36.5 and 37.0 %. By means of a column permeameter and et al.'s [2000] electrokinetic measurement system, we measured a hydraulic permeability between 20 and 22 10 -12 m 2 . We measured the electrical conductivity of the saturated sand σ r for different KCl solutions with conductivity σ ranging between 2.5 and 100 mS m -1 . The electrical formation factor F and the surface conductivity σ S , defined by σ r = σ / F + σ S , were equal to 4.70 and 0.6 mS m -1 , respectively.

Jouniaux

The NE34 (SIFRACO) sand from Nemours, France, is 100 % pure quartz. We did not observe any traces of clays or micas. The grain diameter ranges between 75 and 320 µm, with a mean of 198 µm. The porosity is about 36 % and the permeability 17.5 10 -12 m 2 . The formation factor was found equal to 4.90 (with σ between 50 and 500 mS m -1 ) and there was no detectable surface conductivity, which reflects the absence of clays or micas.

Figure 2 shows the measured fluid conductivity perturbations ∆σ(t) = σ(t) -σ 0 , where σ 0 is the conductivity of the primary fluid, for both the upstream (solid circles obtained from ∆σ(t) = (σ p -σ 0 ) exp (-qt / V up ), an equation derived from mass conservation [START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF], where σ p denotes the upstream fluid conductivity just after the KCl injection, q the flowrate (m 3 s -1 ), t the elapsed time after injection (s), and V up the water volume in the upstream reservoir (m 3 ). In the downstream reservoir, ∆σ started to increase only after the KCl front arrived at the downstream reservoir (at t = t Y ), reached a maximum, and then slowly decreased (Figure 2). The maximum of the downstream ∆σ for the S4M experiment (Figure 2a) occurs at a later normalized time and is 2.5 times smaller than for NE34 (Figure 2b). The significance of this difference will be discussed in section 5.

The results of one KCl front experiment in S4M sand are shown in Figure 3a. The primary fluid was deionized water with a small NaCl content (σ 0 = 8.5 mS m -1 , σ p = 112 mS m -1 and σ p -σ 0 = 103.5 mS m -1 ). The flow rate q was equal to 46 cm 3 min -1 (v = 4.36 10 -5 m s -1 ). Figure 3a shows the EPD curves versus normalized time (here, t Y = 118.4 minutes). We observe a simultaneous and identical drop in all the EPDs at t = 0 to a minimum value around -7 mV followed by a slow increase. Afterwards, the EPDs successively jumped back to zero (and over zero) at sequential times in good agreement with the theoretical arrival times of the front at the electrodes (respectively 0.16, 0.39, 0.61 and 0.84 in normalized time, indicated by the vertical lines in Figure 3). The jumps were followed by broad and large (up to 5 mV) overshoots, whose size increased significantly with distance from the upstream reservoir, and ending with a slow return to a stabilized value around -1 mV. Other S4M experiments were performed using tap water (σ 0 ≈ 45 mS m -1 , σ p ≈ 145 mS m -1 ). The EPDs exhibited a very similar shape, although the amplitude of the negative part was smaller (around -5 mV).

For the pure quartz sand experiment shown in Figure 3b, tap water served as primary fluid, with σ 0 = 46 mS m -1 , σ p = 143 mS m -1 and q = 34.6 cm 3 min -1 (thus v = 3.4 10 -5 m s -1 , t Y = 151.8 min). The initial drop of the EPDs (Figure 3b) had a magnitude more than 3 times smaller than that in the S4M sand. The decrease rate was rapid at first and then slowed down with time. As before, the EPDs showed a sequence of very sharp increases corresponding to the theoretical arrival times of the concentration front at each particular electrode, followed by a small, very narrow overshoot and, finally, a rapid return to zero. Another NE34 experiment (q ≈ 42 cm 3 min -1 , same σ 0 and σ p ) gave similar results.

DISCUSSION

Maineult et al.

[2004] devised a conservative model (i.e., without chemical reactions) to compute the downstream conductivity variation ∆σ. The results of this model are represented by the dotted lines in Figure 2. Experimental data and model agree quite well for the pure quartz sand only (Figure 2b). For S4M sand, the model overestimates the maximum by a factor 2 and underestimates its time of occurrence. This observed deficit in conductivity suggests that potassium chloride interacted with some of the minor phases contained in S4M sand but not present in NE34 sand. Experiments on S4M using NaCl [START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF] exhibited a similar but much smaller effect (the predicted maximum of ∆σ in the downstream reservoir was less 20 % higher than observed); moreover, in the NaCl case, the SP signals did not exhibit positive humps as in Figure 3a. Therefore, we rule out Cl -as main reacting species and conclude that the chemical processes had to involve the K + ions. [START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF][START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF] constructed a conservative model for the SP response to the advective transport of concentration fronts (note that the fluid concentration in the sand-body and its evolution with time must first be determined [for details see [START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF][START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF]). This model was validated by comparison with the results of NaCl advection experiments in S4M sand. As briefly explained in section 2, the electrical potential difference ∆U(y e ,t) between an electrode located at y e and the reference electrode is the sum of two terms:

( ) ( ) ( ) ( )
∆ =∆ + ∆ -∆ ⎡ ⎤ ⎣ ⎦ (4)
where ∆U j (y e ,t) denotes the junction potential contribution, which is obtained by integrating equation ( 2) between 0 and y e . The bracketed term corresponds to the electrokinetic contribution and is obtained by integration of equation ( 1):

( ) ( ) ( ) ( ) ( ) * * 0 , , 0

∫

(5) (Note that the incorrectly typed equation ( 11) in [START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF] should be identical to this last equation). Numerical computation of the integral in equation ( 5) requires knowledge of the relation L * (σ(y,t)), which we experimentally estimated on S4M and NE34 samples using the device and procedure described in [START_REF] Jouniaux | Streaming potential in volcanic rocks from Mount Pelée[END_REF]. For each sand, we used as circulating solutions the same primary fluid as in the advection experiments described above (i.e., tap water for NE34 and deionized + 30 mg L -1 NaCl for S4M) to which 0, 250 or 500 mg L -1 KCl were added. Prior to the measurements, we circulated these solutions in a closed-circuit for a sufficiently long time to reach chemical equilibrium. In accordance with equation ( 1), which predicts that L * is inversely proportional to σ, we determined that L * = -1.7 10 -4 σ -1 in the range 50 to 150 mS m -1 for NE34 and L * = -2.8 10 -4 σ -1 in the range 5 to 100 mS m -1 for S4M. Then, after solving the conservative transport problem for KCl using the method of [START_REF] Maineult | Electrical response of flow, diffusion and advection in a laboratory sand-box[END_REF][START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF], we computed the EPDs from equation (4) for the NE34 sand (Figure 3c). A perfect match is not obtained. The model predicts sharp, angular EPD curves at t = 0 whereas the observed ones are rounded. Also, the small, narrow overshoots are not reproduced in the theoretical EPDs. However, the model successfully predicts the moderate amplitude of the signals, suggesting that, to first order, the EPDs are dominated by the junction potential (the electrokinetics contribution, not shown here, represents less than 2 % of the total potential). Note that the initial rounding of the EPDs, not predicted by our model, was also observed in all our previous experiments when the reference electrode was placed in the upstream reservoir, irrespective of the electrolyte used [e.g., [START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF]. The general occurrence of this effect suggests that it is more likely caused by insufficient precision of the concentration modeling rather than by the SP model itself. In the aforementioned previous experiments, we also always observed narrow peaks accompanying the passage of the concentration front (just before or just after) but their magnitude (0.3 to 0.6 mV) was much smaller than the overshoots of Figure 3a.

Application of our conservative SP model to S4M (using the appropriate value for L * ) produces a response similar to those shown in Figure 3c and thus disagreeing completely with the observations. The fact that the model is somewhat successful in the case of pure quartz sand NE34 but fails totally in S4M suggests that the minor mineralogical phases present in S4M (but absent in NE34) must be responsible for this effect.

Based on the two observations above (i.e., deficit of downstream conductivity and failure of the SP model for the S4M sand) we propose that chemical reactions with the minor minerals in S4M sand likely occurred, leading to trapping of K + during the passage of the concentration front and their release during passage of the tail. Prior examination of the virgin S4M sand using a binocular magnifying glass showed that the mica grains were essentially biotite, and, to a much smaller degree, muscovite. A new examination, carried out after the NaCl transport experiments of Maineult et al.'s [2004Maineult et al.'s [ , 2005] ] were performed (but before the KCl ones), revealed that the biotite grains had been altered. It is well known during weathering alteration, biotite is chloritized by potassium release [e.g., [START_REF] Mitchell | Potassium and argon loss patterns in weathered micas: implications for detrital mineral studies, with particular reference to the triassic palaeogeography of the british isles[END_REF]1984, and references herein] and that this reaction is quite rapid [START_REF] Craw | Oxidation and microprobe-induced potassium mobility in iron-bearing phyllosilicates from the Otago schists[END_REF]Pal, 1985;[START_REF] Malmström | The dissolution of biotite and chlorite at 25°C in the near-neutral pH region[END_REF]. During a biotite dissolution experiment in neutral pH conditions, [START_REF] Malmström | The dissolution of biotite and chlorite at 25°C in the near-neutral pH region[END_REF] observed that the K + concentration in the fluid remained constant at an equilibrium value (1 mmol L -1 ). The stability of the K + concentration was explained by a dynamic equilibration between trapping and releasing of K + . In our experiment, the KCl concentration at the front was equal to about 7 mmol L - 1 , i.e. greater than the equilibrium concentration mentioned above. Thus, we imagine that, in experiments reported here, K + ions were trapped by chloritized biotites as long as the KCl concentration remained over this critical concentration (around 1 mmol L -1 ). The cationic conductivity σ + in equation ( 1) therefore decreases and the (negative) potential gradient increases in absolute value, corresponding to the observed enhanced amplitude of the initial drop (compared to the response in pure quartz sand). K + trapping also explains the deficit of conductivity in the downstream reservoir. Then, as the front concentration passes below the equilibrium concentration, the trapped potassium ions should be released, leading to an increase in σ + and to the consequent change in polarity, i.e., to the big positive humps observed in the SP signals. This releasing may explain the fact that the maximum of the downstream conductivity was delayed in time (Figure 3a).

This study experimentally demonstrates that chemically active, mineral phases present in rocks and soils even in small amount can strongly affect the SP response to transport of ionic fronts, although we have not yet determined the precise chemical mechanisms involved in the particular case studied here. The retardation effect illustrated here predominantly affects the junction potential term in the SP signals. In other situations, it has been shown that chemical reactions can modify the electrokinetic term as well. For example, the electrokinetic coupling coefficient was significantly reduced in absolute value, and even changed in sign, during calcite precipitation in sand [START_REF] Guichet | Modification of streaming potential by precipitation of calcite in a sand-water system: laboratory measurements in the pH range from 4 to 6[END_REF]. Chemical modifications of both fluid and matrix, such as those induced by dissolution/re-precipitation at high temperature, can also strongly affect the SP response [START_REF] Bernabé | Laboratory measurements of electrical potential in rocks during high-temperature water flow and chemical reactions[END_REF]. The electrokinetic behavior of rocks containing different mineral phases (but without any reactions) can also be quite complicated [START_REF] Hase | ζ potential measurements of volcanic rocks from Aso caldera[END_REF].

Despite their importance for field data interpretation (it seems highly probable that chemically active minerals will be present in soils and rocks at any site of interest), the effects of the mineralogical composition of rock and of the chemical reactions between minerals and fluid on the self-potential still remain poorly documented. Many studies focused on the electrokinetic potential of homogeneous media (i.e., nearly pure quartz sand or pure clay powder), and few experiments have been performed to date to study the junction potential in porous media. There is presently not enough information available to accurately predict the SP response in strongly heterogeneous soils and rocks. The incremental approach, i.e., starting from an already well-known system and modifying it gradually, may help to gain understanding on this complex topic.

  ) and downstream (open circles) reservoirs during two experiments using the S4M and NE34 sands. Note that to compare the results of the different experiments (the sand permeability and, therefore, the flowrate varied between runs owing to slight differences in sand compaction), we normalized the time by dividing it by the time needed by the front to reach the downstream reservoir (neglecting the hydrodynamic dispersion) t Y = vY, where Y = 31 cm, and ν is the mean fluid velocity. In both experiments, the values measured in the upstream reservoir are in good agreement with theoretical values (continuous lines)

Figure 1 .

 1 Figure 1. Experimental set-up. Flow is generated by tilting the sandbox while

Figure 2 .

 2 Figure 2. Observed (dots) and predicted (lines) variations of the fluid electrical

Figure 3 .

 3 Figure 3. Measured EPDs versus normalized time for the S4M (a) and NE34 (b)

K -1 ), A Avogadro's number (≈6.022 10

mol -1 ) and e the absolute electron charge


Andrés Galavís

David González

-Elias Cueto
email: ecueto@unizar.es

Francisco Chinesta

Manuel Doblaré

A natural element updated Lagrangian approach for modelling fluid structure interactions

Keywords: fluid structure interaction, meshless methods, natural element method, -shapes. MOTS-CLÉS : fluide structure interaction, méthode sans maillage, méthode des éléments naturels, formes-1

In this paper we present a novel methodology for the numerical simulation of fluid structure interactions in the presence of free surfaces. It is based on the use of the Natural Element Method (NEM) in an updated Lagrangian framework, together with the integration of the Navier-Stokes equations by employing a Galerkin-characteristics formulation. Tracking of the free-surface is made by employing shape constructors, in particularshapes. A theoretical description of the method is made and also some examples of the performance of the technique are included.

RÉSUMÉ. Dans cet article on présente une nouvelle méthodologie pour la simulation numérique des interactions fluide-structure en présence de surfaces libres. Elle se base sur l'emploi de la méthode des éléments naturels (MEN) dans un cadre Lagrangien actualisé, utilisant en même temps l'intégration des équations de Navier-Stokes à travers la méthode des caractéristiques-Galerkin. Des constructeurs de forme sont employés pour le suivi de la surface libre, particulièrement des formes-. On présente une description théorique de la méthode, ainsi que quelques exemples du fonctionnement de la technique.

Introduction

The fact that meshless methods [START_REF] Belytschko | Element-Free Galerkin Methods[END_REF][START_REF] Liu | Reproducing kernel particle methods[END_REF] do not suffer of mesh distortion opened a renewed interest in the last decade in Lagrangian formulations for some problems, being free-surface flows a typical example. Thus, it is possible to employ an updated Lagrangian strategy for the fluid domain, while employing a total or updated Lagrangian strategy for the solid. This approach is very convenient for some classes of problems, especially those involving drastic changes in the fluid domain geometry. Both domains are then formulated in similar frameworks and the coupling between them becomes more direct than in ALE formulations (see [START_REF] Donea | Arbitrary Lagrangian-Eulerian finite element methods[END_REF] or [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF]).

In this paper we describe mainly the fluid flow formulation proposed in the context of an updated Lagrangian strategy. We employ the -shape-based Natural Element Method ( -NEM) [START_REF] Cueto | Modeling three-dimensional piece-wise homogeneous domains using the -shape based Natural Element Method[END_REF]2003) to this end. At present the solid is assumed rigid, being prescribed its kinematics. This formulation posses some advantages, that include an exact interpolation along the boundary (Cueto et al., 2001), that allows for a standard, FE-like, treatment of the fluid-solid interface conditions. We firstly describe the bases of the -NEM and then introduce the proposed numerical scheme for the integration of the Navier-Stokes equations. Finally, we include some examples that demonstrate the accuracy of the proposed scheme and also prove the potential of the technique.

The natural element method

Standard formulation

The NEM [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF]Cueto et al., 2003) is a Galerkin procedure based on the natural neighbor interpolation scheme, which in turn relies on the concepts of Voronoi diagrams and Delaunay triangulations (see Figure 1), to build Galerkin trial and test functions. These are defined as the Natural Neighbor coordinates of the point under consideration, that is, with respect to Figure 2, the value in the point x of the shape function associated with the node 1 is [START_REF] Sibson | A Vector Identity for the Dirichlet Tesselation[END_REF][START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF]. In addition, the NEM has other interesting properties such as linear consistency and smooth shape functions (C 1 everywhere except of the nodes). These functions are dependent on the position and density of nodes, leading to standard FE constant strain triangle shape functions, bilinear shape functions or rational quartic functions in different situations (see Figure 3 for a typical shape function). These properties permit an exact reproduction of linear displacement fields on the boundary of convex domains.

abcd abfe A A x ) ( 1 [1]

Figure 2. Definition of natural neighbour coordinates

-shape formulation

A slight modification of the way in which the Natural Neighbour interpolant is built was proposed to achieve linear interpolation also over non-convex boundaries (Cueto et al., 2001). This modification was based on the concept of -shapes. These are a generalization of the concept of the convex hull of a cloud of points and are widely used in the field of scientific visualization and computational geometry to give a shape to a set of points. Alpha-shapes give shape to a cloud of points and are widely used in Computational Geometry despite having been developed quite recently. They were first introduced in two-dimensions by Edelsbrunner in 1983, and not generalised in three-dimensions until (Edelsbrunner and Mücke, 1994). Anshape is a generalisation of the convex hull of a cloud of points. It is a polytope that is not necessarily convex and that can be triangulated by a subset of the Delaunay triangulation, thereby maintaining the empty circumcircle criterion.

In what follows, we introduce the formal definition of a complete family ofshapes for a given set of points N, as in Edelsbrunner and Mücke (1994). Let N be a finite set of points in 3 and a real number with 0

. A k-simplex T with 0 k 3, is defined as the convex hull of a subset T N of size T =k+1. Let b be an -ball, that is, an open ball of radius . A k-simplex T for 0 k 2 is said to beexposed if there exists an empty -ball b with T = b N, where indicates the boundary of the ball or, more properly, the sphere or plane bounding b. That is, a ksimplex is -exposed if an -ball whose boundary passes through its defining points contains no other point of the set N. In this way, we can define a family of sets F k, as the sets of -exposed k-simplices for the given set N, fixed and 0 k 2.

Based on these concepts, the -shape of N, S , is defined as the polytope whose boundary consists of the triangles in F 2, , the edges in F 1, and the points or vertices in F 0, . As the value decreases, the -shape shrinks by the progressive development of cavities or holes. For this to occur, one or more -balls can occupy the interior of a simplex. The value clearly gives an intuitive measure of the maximum curvature in a region of the domain. The -shape concept is also a generalisation of the convex hull since the -shape for value = 0 is identical to the initial set of points, i.e., S0 = N, and the -shape for sufficiently high values of is the convex hull of the given set.

An example of some -shapes of the complete family for a given set of points distributed over the geometry of a human jaw can be seen in Figure 4. It has been demonstrated (Cueto et al., 2001) how the construction of the interpolant over an appropriate -shape of the domain gives rise to an exact imposition of essential boundary conditions over any kind of domain (convex or not.) In addition, it enables us to track the flow front position accurately.

A natural neighbour updated Lagrangian strategy for the fluid domain

In this section we review the time integration scheme developed in (Gonzalez et al., 2006), that will be applied in the integration of the fluid flow equations. It is based on a Galerkin-characteristics formulation of the Navier-Stokes equations.

Governing equations

We consider here the problem of Fluid Dynamics at moderate Reynolds number. Thus, the governing equations can be set as follows. Consider a fluid in a region of the space 2 or 3 . The fluid flow is governed by the following momentum and mass balance equations:

b v v v ) ) • ( ( ,t in t , 0 [2] 0 • v in t , 0 [3]
where v represents the fluid velocity, the stress tensor, represents fluid density and b the volumetric forces acting on the fluid. The constitutive equation for a newtonian fluid is given by:

D I I 2 p p [4]
where D is the strain rate tensor, p the pressure and the dynamic viscosity of the fluid. To solve the problem we must prescribe an initial state as well as boundary conditions, as usual.

Time discretisation

The motion equations can be grouped to

v v v v b ) • ( • t dt d [5] 0 • v [6] D I 2 p [7]
The weak form of the problem associated to Equations The second term in the right-hand side of Equation [8] represents the inertia effects. Time discretization of this term represents the discretization of the material derivative along the nodal trajectories, which are precisely the characteristic lines related to the advection operator. Thus, assuming known the flow kinematics at time t n t n ) 1 ( 1, we proceed as follows:

d t d dt d n n * 1 * ) ( ) ( v X v x v v v [9]
where X represents the position at time t n-1 occupied by the particle located at position x at present time t n , i.e.:

t n ) ( 1 X v X x [10]
So we arrive to:

d t d d t d p d n * 1 * * * * : : 2 v v bv vv D I D D [11]
where we have dropped the superscript in all the variables corresponding to the current time step.

Algorithmical issues

The most difficult term in Equation [11] is the second term of the right-hand side. The numerical integration of this term depends on the quadrature scheme employed. If we employ traditional Gauss-based quadratures on the Delaunay triangles, it will be necessary to find the position at time t n-1 of the point occupying at time t n the position of the integration point k :

k k k k n n t d t ) ( ) ( * 1 * 1 v v v v [12]
where k represent the weight associated to integration point k , and k corresponds to the position occupied at time t n-1 by the quadrature point k . If we employ some type of nodal integration, as in [START_REF] Jiun-Shyan | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF]Gonzalez et al., 2004a), this procedure becomes straightforward, with the only need to store nodal velocities at time step t n-1 .

We discuss here the procedure to follow when employing Gauss quadratures on the Delaunay triangles. We proceed iteratively. Denoting by i the current iteration, ( 1 i ), we apply

t i k n i k k ) ( 1 1 X v X x with 0 k k X x [13] until 1 i k i k X X
within a prescribed tolerance.

We have assumed that the number of natural neighbours of a given integration point does not change during a time step, thus needing the storage of nodal velocities at time t-1 only. It can occur that some of the nodes neighbouring the integration point at time t were not actually its neighbours at time t-1, but this does not constitute a problem, since the number of natural neighbours of a point is usually high (much bigger than three), so the quality of the interpolation is thus guaranteed. In fact, this procedure has shown to converge at a high speed, with no more than 3 iterations, at least for reasonable time steps.

Numerical examples

Broken dam problem

The broken dam problem is classic when testing the performance of integration methods for free surface flows. We consider a rectangular column of water, initially retained by a door that is instantaneously removed at time t=0 (see Figure 5). When the door is removed, water flows under the action of gravity, considered as 9.81 m/s 2 . Density of water is 10 3 kg/m 3 , and a viscosity of 0.1 Pa•s was assumed as in other numerical simulations performed using different numerical strategies (see [START_REF] Duchemin | Jet formation in bubbles bursting at a free surface[END_REF] and references therein, for instance). The discrete model was composed of 3 364 nodes. No remeshing, addition or deletion of nodes was performed throughout the computation.

Figure 5. Geometry of the broken dam problem

Figure 6 shows a comparison between numerical results and experimental ones, obtained from the literature [START_REF] Martin | Part IV. an experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF]. As can be noticed, an excellent agreement was found between experimental and numerical results, despite the distortion of the triangulation. In Figure 7 the error in mass conservation is depicted, which remained always below 3%. The influence of the relationship between the parameter and the nodal parameter h on this error was deeply analyzed in [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF]. Stick boundary conditions were assumed on the reservoir walls, being the upper water surface a free boundary which evolves slightly during the simulation as noticed in Figure 8c. Time increment was set to 0.005s. being the fluid viscosity of 0.01 Pa.s.

The ability of the proposed method for describing flows in the framework of updated Lagrangian description is then fully proved.

Water mill partially submerged

The proposed method seems particularly well adapted for dealing with freesurface flows. If the sail is only partially submerged, then large-amplitude waves are expected, justifying the interest of the present simulation. For this purpose, we consider the same geometry as in the previous example, but maintaining the sail only partially submerged, as shown in Figure 9.

Figure 9. Initial geometry of the mill problem for the partially submerged configuration

Material parameters were chosen as in the previous example. In this case the time step was set to 0.03s. This test can be found in other references, see, for instance [START_REF] Idelsohn | The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF]. Note the appearance of a large amplitude wave on the free surface of the liquid. The geometry of the fluid and the eventual generation of drop and jets can be accurately described by the -shapes. A deep study on this topic has been recently presented in [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF], Gonzalez et al., 2006).

Conclusions

This paper proposes a Galerkin-characteristics updated-Lagrangian fluid flow formulation for simulation of fluid structure interaction problems. The fields approximation is based on the use of the Natural Element Method which makes it possible to work with the same cloud of nodes which moves with the material velocity, avoiding remeshing stages. The use of Lagrangian descriptions in both the solid and fluid domains greatly simplifies the formulation and numerical resolution of fluid structure interaction problems, especially those involving free-surfaces.

The application of the proposed scheme to real FSI problems (i.e., those in which the movement of the solid is coupled with the fluid one) is currently the aim of our research.

Figure 1 .

 1 Figure 1. Delaunay triangulation and Voronoi diagram of a set of points

Figure 3 .

 3 Figure 3. Typical Sibson shape function (courtesy N. Sukumar)

Figure 4 .

 4 Figure 4. Five elements of the complete family of -shapes of a cloud of points distributed over the geometry of a human jaw. Increasing values of from 0 (top line, left) to infinity (bottom line, right)

Figure 6 .Figure 7 .

 67 Figure 6. Front position (in non-dimensional form) in time. Numerical results vs. experimental ones

Figure 10 .

 10 Figure 10. Six snapshots of the generation of a wave during the rotation of the mill

Acknowledgement

The support from the Spanish Ministry of Education and Science, through project DPI-2005-08727-C02-01 is gratefully acknowledged.