Claude Lefèvre
email: clefevre@ulb.ac.be

Stéphane Loisel
email: stephane.loisel@univ-lyon1.fr

Claude Lef Èvre

ON FINITE-TIME RUIN PROBABILITIES FOR CLASSICAL RISK MODELS

Keywords: ruin probability, finite and infinite horizon, compound binomial model, compound Poisson model, ballot theorem, pseudo-distributions, Solvency II, Value-at-Risk

 for the probability of (non-)ruin within finite time. First, a standard method based on the ballot theorem and an argument of Seal-type provides an initial (known) formula for that probability. Then, a concept of pseudo-distributions for the cumulated claim amounts, combined with some simple implications of the ballot theorem, leads to the desired formula. Two expressions for the (non-)ruin probability over an infinite horizon are also deduced as corollaries. Finally, an illustration within the framework of Solvency II is briefly presented.

Introduction

The two classical models in risk theory are the compound binomial model, which is a discrete-time process where the claim amounts are usually assumed to be integer-valued random variables, and the compound Poisson model, which is a continuous-time analogue of this process where the claim amounts are generally assumed to have an absolutely continuous distribution. Although the continuous-time version is quite popular, a discretetime version can be appropriate for some applications and provides a better intuitive understanding.

The present paper deals with the problem of evaluating, for both risk models, the probability of (non-)ruin within finite time. This problem has received much attention in the literature, and different algorithms have been proposed, with their own advantages and drawbacks; see, e.g., the analysis by Dickson [START_REF] Dickson | On numerical evaluation of finite time survival probabilities[END_REF] and the references therein.

Recently, Picard and Lefèvre [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] derived an elegant explicit formula, called P.L. formula below, for the finite-time non-ruin probability in a compound Poisson model where the claim amounts are integer-valued. Such a case is important because in practice a discretization of the claim amounts is often required for numerical calculations (e.g., de Vylder and Goovaerts [START_REF] De | Recursive calculation of finite-time ruin probabilities[END_REF]). The importance of the P.L. formula has been pointed out by De Vylder [START_REF] De | La formule de Picard et Lefèvre pour la probabilité de ruine en temps fini[END_REF], [START_REF] De | Numerical finite-time ruin probabilities by the Picard-Lefèvre formula[END_REF] and Ignatov et al. [START_REF] Ignatov | An improved finite-time ruin probability formula and its Mathematica implementation[END_REF]. A refined discussion of its merits in comparison with other available formulas is given in Rullière and Loisel [START_REF] Rullière | Another look at the Picard-Lefèvre formula for finite-time ruin probabilities[END_REF]. Besides, De Vylder and Goovaerts [START_REF] De | Explicit finite-time and infinite-time ruin probabilities in the continuous case[END_REF] proved that a similar formula holds too for a compound Poisson model with continuous claim amounts; see also Ignatov and Kaishev [START_REF] Ignatov | A finite-time ruin probability formula for continuous claim severities[END_REF].

Several methods are possible to establish the P.L. formula. In Picard and Lefèvre [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] and De Vylder and Goovaerts [START_REF] De | Explicit finite-time and infinite-time ruin probabilities in the continuous case[END_REF], the first step consists in stipulating for a total claim amount until some time t and conditioning on the last claim instant and amount before t. This yields an integral equation which is then solved by the former authors using the generalized Appell polynomials (in an extended framework), and by the latter authors using the convolution products and the Laplace transforms. On another hand, Rullière and Loisel [START_REF] Rullière | Another look at the Picard-Lefèvre formula for finite-time ruin probabilities[END_REF] showed that the P.L. formula can be linked to the well-known ballot theorem and a Seal-type formula (see, e.g., Seal [START_REF] Seal | The Stochastic Theory of a Risk Business[END_REF] and Gerber [START_REF] Gerber | An Introduction to Mathematical Risk Theory[END_REF]).

The primary aim of the present paper is to establish that the formula derived by Picard and Lefèvre [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] can be generalized to the classical compound binomial and compound Poisson risk models. As in its original version, this extended formula will be given explicitly in terms of the pseudo-distributions of some cumulated claim amounts. The method of proof will enable us to provide a unified approach to the ruin problem over a finite horizon. To begin with, we obtain a first formula for the finite-time probability of non-ruin by applying a standard argument, based on the ballot theorem when the initial reserves are null and on a Seal-type formula when the reserves are positive. The result here is not new but this will put us in a position to develop an alternative argument. Succinctly, we first introduce a concept of pseudo-distributions for the cumulated claim amounts, and we show how to compute numerically these pseudo-laws using simple recursions. This mathematical tool, combined with some simple implications of the ballot theorem, will then enable us to derive the desired formula of P.L.-type. In addition, we will also deduce from the Seal and P.L.-types formulas two expressions for the (non-)ruin probability over an infinite horizon. Finally, a brief illustration is presented that points out how finite-time ruin probabilities could be of interest within the framework of Solvency II.

Compound binomial risk model

This Section is concerned with a discrete-time risk model equivalent to the compound binomial model, but given under a less traditional formulation where more than one claim can arise per time period. We mention that such a model has been investigated by Gerber [START_REF] Gerber | Mathematical fun with the compound binomial process[END_REF], De Vylder and Goovaerts [START_REF] De | Recursive calculation of finite-time ruin probabilities[END_REF], Shiu [START_REF] Shiu | The probability of eventual ruin in the compound binomial model[END_REF], Willmot [START_REF] Willmot | Ruin probabilities in the compound binomial model[END_REF], Dickson [START_REF] Dickson | On numerical evaluation of finite time survival probabilities[END_REF], De Vylder and Marceau [START_REF] De | Classical numerical ruin probabilities[END_REF], Cheng et al. [START_REF] Cheng | Discounted probabilities and ruin theory in the compound binomial model[END_REF], Li and Garrrido [START_REF] Li | On the time value of ruin in the discrete time risk model[END_REF] and Picard et al. [START_REF] Picard | Problèmes de ruine en théorie du risque à temps discret avec horizon fini[END_REF].

Let us consider an insurance company that evaluates its reserves at periodic times fixed in advance. Claims arise at any date but they can be registered only at the end of the period during which they occur. The successive claim amounts per period, X i for i = 1, 2, . . ., are non-negative integer-valued independent identically distributed (i.i.d.) random variables; let {p n , n ∈ IN} be their common probability mass function (IN ≡ {0, 1, 2, . . .}). At the beginning of each period, the company receives a constant premium, Figure 1: Cumulated premium and claim amounts. assumed to be equal to 1, as a prepayment to cover the risk during that period (see the remark below). The initial reserves are of amount u ∈ IN.

Take the length of a period as the time unit, and denote by t ∈ IN the new time scale. During any time interval (0, t], the cumulated claim amount is given by S t = X 1 +. . .+X t , and the total premium income is equal to t. Thus, the reserves R t of the company at time t are given by

R t = u + t -S t , t ≥ 1. (2.1)
Note that at time t + 0, the company will receive an additional premium of 1 but this is a prepayment for the next period (t, t + 1]. A typical situation is represented in Figure 1.

Ruin will occur at the first time T u when the reserves become negative or null. In other words, T u ≥ 1 necessarily and for t ≥ 1,

T u ≥ t + 1 means S i < u + i for 1 ≤ i ≤ t. (2.2)
This definition of ruin is that adopted by Gerber [START_REF] Gerber | Mathematical fun with the compound binomial process[END_REF] rather than Shiu [START_REF] Shiu | The probability of eventual ruin in the compound binomial model[END_REF] where ruin occurs when the reserves become negative. The exact distribution of T u is the topic of investigation hereafter. Let us mention that if the premium per time unit is equal to a positive constant c = 1, so that the reserves are R t = u + ct -S t , the ruin problem is different from the actual one. Indeed, the monetary unit might be adapted by dividing by c, but this would give claim amounts that are no longer integer-valued. Nevertheless, the ruin probability for this case can be evaluated by following the method developed in Picard et al. [START_REF] Picard | Problèmes de ruine en théorie du risque à temps discret avec horizon fini[END_REF].

A standard approach

Basic case u = 0. Ballot-type problems are among the oldest questions discussed in probability. Their applications are numerous, in risk theory as here and in many other fields such as queueing, reliability, sequential analysis and random graphs. Much on this can be found in the book by Takács [START_REF] Takács | Combinatorial Methods in the Theory of Stochastic Processes[END_REF]. A well-known generalization of the original ballot theorem due to Takács [START_REF] Takács | A generalization of the ballot problem and its application in the theory of queues[END_REF] is as follows.

Lemma 2.1 (Ballot theorem, discrete case). Let X i , i ≥ 1, be a sequence of non-negative integer-valued i.i.d. random variables. Consider the associated partial sums process S i = X 1 + . . . + X i , i ≥ 1. Then, for any integers t ≥ 1 and 0 ≤ n ≤ t -1, P (S i < i for 1 ≤ i < t, and

S t = n) = t -n t p * t n , (2.3)
where {p * t n , n ∈ IN} denotes the probability mass function of S t (i.e. the t-th convolution of the law of X 1 ).

Graphically, the left-hand side of (2.3) represents the probability that the partial sums process reaches the level n at time t without meeting or crossing earlier the diagonal line of slope 1 (which is an upper boundary); see Figure 2.

Note that formula (2.3) remains true for exchangeable sequences. Many proofs of this result have been proposed, including methods based on a combinatorial argument (e.g., Takács [START_REF] Takács | A generalization of the ballot problem and its application in the theory of queues[END_REF]), a stopping time theorem for backward martingales (e.g., Grimmett and Stirzaker [START_REF] Grimmett | Probability and Stochastic Processes[END_REF]) or a queueing theoretic construction (Konstantopoulos [20]). To get a self-contained presentation, we choose to give below an elementary proof in the i.i.d. case, by applying a simple trick that will be useful too later.

Proof of (2.3). Let us make a rotation of 180 • in Figure 2, and consider the partial sums trajectories in reversed time. Since the increments of the process are i.i.d., each sample path keeps the same probability after rotation. Note also that the reversed partial sums process starts at level t -n and is backward shifted in the sense that the first jump X 1 occurs at time 0, the second jump X 2 at time 1, and so on; see Figure 3. Therefore, (2.3) is equivalent to the following identity: for t ≥ 1 and 0 ≤ n ≤ t -1, P (t -n + S i > i for 1 ≤ i < t, and t -n + S t = t) = t -n t p * t n .

(2.4)

For m ≥ 1, let τ m be the first-meeting time of the process m + S i , i ≥ 1, with the diagonal line (which is a lower boundary). We observe that the probability in the lefthand side of (2.4) corresponds to P (τ t-n = t). Obviously, τ m ≥ m and P (τ m = m) = p * m 0 , m ≥ 1. For n ≥ 1, τ m = m + n means that firstly, S m takes an arbitrary integer value j say, between 1 and n, and then, independently, a similar shifted partial sums process starting at level j meets the diagonal line at time n (Figure 3); this gives

P (τ m = m + n) = n j=1 p * m j P (τ j = n).
(2.5)

We point out that (2.5) provides a recursive formula to calculate the probabilities

P (τ m = m + n) for m, n ≥ 1.
On another hand, let us consider the partial sums S m , m ≥ 1. The assumption of i.i.d. increments implies that, for n ≥ 1,

E(S m | S m+n = n) = 1 p * (m+n) n n j=1 j p * m j p * n n-j , (2.6)
Figure 2: Non-crossing of {S i } through the diagonal line y = x. as well as

E(S m | S m+n = n) = m n m + n . (2.7)
Combining (2.6) and (2.7) then yields the identity

m m + n p * (m+n) n = n j=1 p * m j j n p * n n-j . (2.8)
Let us now define the quantities

θ m (m + n) ≡ [m/(m + n)] p * (m+n) n , for m ≥ 1, n ∈ IN. Thus, θ m (m) = p * m 0 and for n ≥ 1, (2.8) gives θ m (m + n) = n j=1 p * m j θ j (n).
(2.9)

Here too, (2.9) allows us to determine recursively the quantities θ m (m + n). In fact, the recursive formulas (2.5) and (2.9) are identical, so that

P (τ m = m + n) = θ m (m + n).
Formula (2.3) then follows by putting m = t -n. ⋄

Let us turn to the ruin time T 0 for the risk model with u = 0. Of course, T 0 = 1 when X 1 ≥ 1, and for t ≥ 1, we can write that

P (T 0 ≥ t + 1) = t-1 n=0 P (S i < i for 1 ≤ i ≤ t, and S t = n).
(2.10)

Applying the ballot formula (2.3) to the right-hand side of (2.10) then leads to the following result (e.g., Takács [START_REF] Takács | Combinatorial Methods in the Theory of Stochastic Processes[END_REF]).

Proposition 2.2 (Takács-type formula). P (T 0 = 1) = P (X 1 ≥ 1), and

P (T 0 ≥ t + 1) = 1 t t-1 n=0 (t -n) p * t n , t ≥ 1. (2.11)
We notice that the probability of non-ruin before time t is given explicitly in terms of the distribution of S t evaluated at points n = 0, . . . , t -1.

Passage to u ≥ 0. Looking at Figure 2, let us now suppose that the upper boundary for the partial sums process S i , i ≥ 1, is a straight line of slope 1 beginning at level u ∈ IN instead of 0.

Lemma 2.3 For any integers

u ≥ 0, t ≥ 1 and u ≤ n ≤ u + t -1, P (S i < u + i for 1 ≤ i < t, and S t = n) = p * t n - n j=u+1 t + u -n t + u -j p * (t+u-j) n-j p * (j-u) j .(2.12)
Proof. Operating a rotation of 180 • as before, we observe that

P (S i < u + i for 1 ≤ i < t, and S t = n) = P (τ t+u-n ≥ t and S t = n), (2.13)
where τ t+u-n denotes again the first-crossing time of the process t + u -n + S i , i ≥ 1, with the (lower) diagonal line. Now, for any m ≥ 1, P (m + S t = m + n) = p * t n . On another hand, the process m + S i , 1 ≤ i < t, will either cross or not the diagonal line, and if yes, for the first time at some level k say, between m and t -1. Therefore, we can write that

P (m + S t = m + n) = P (τ m ≥ t and S t = n) + t-1 k=m P (τ m = k) P (S t-k = m + n -k). (2.14)
Using formula (2.3), we obtain from (2.14) that

P (τ m ≥ t and S t = n) = p * t n - t-1 k=m m k p * k k-m p * (t-k) m+n-k . (2.15)
Putting m = t + u -n in (2.15) and taking j = t + u -k as a new index of summation then yields formula (2.12). ⋄

Remark. When u = 0, (2.12) and (2.3) hold true and the following identity then follows: for t ≥ 1 and 0

≤ n ≤ t -1, p * (t) n = n j=0 t -n t -j p * (t-j) n-j p * j j . (2.16)
This is precisely a generalization of (2.16) which will be the key result in the analysis of Section 2.3. Formula (2.16) is easily obtained on basis of a graphical representation. Indeed, consider a line of slope 1 passing at point (t, n) (of equation y = n -t + x). The event S t = n implies that the process S i , 1 ≤ i ≤ t, necessarily crosses this line, for the first time at some point (t -j, j) say, with 0 ≤ j ≤ n. Note that the first-crossing time, inf{i : S i = n -t + i}, corresponds to the time τ t-n above. Thus, we have

P (S t = n) = n j=0 P (τ t-n = t -j) p * j j .
which gives (2.16) by virtue of (2.3). ⋄ Let us examine the ruin time T u for the risk model with u ≥ 0. Of course, T u = 1 when X 1 ≥ u + 1, and for t ≥ 1, we get

P (T u ≥ t + 1) = u n=0 P (S t = n) + u+t-1 n=u+1 P (S i < u + i for 1 ≤ i < t, and S t = n). (2.17)
Using (2.12) in the right-hand side of (2.17) and permuting the two sums, we then find formula (2.18) below. A formula of this type was already obtained by Seal [START_REF] Seal | The Stochastic Theory of a Risk Business[END_REF].

Proposition 2.4 (Seal-type formula). P (T u = 1) = P (X 1 ≥ u + 1), and

P (T u ≥ t + 1) = u+t-1 j=0 p * t j - u+t-1 j=u+1 p * (j-u) j u+t-1 n=j t + u -n t + u -j p * (t+u-j) n-j , t ≥ 1. (2.18)
Note that formula (2.18) is given explicitly in terms of the distribution of S 1 , . . . , S t at different points between 0 and u + t -1.

An alternative approach

Pseudo-distributions for S t . To begin with, we are going to point out an underlying polynomial structure in the distribution of S t , t ∈ IN. Without real loss of generality, it is assumed that 0 < p 0 = P (X 1 = 0) < 1. Now, among the t claims with sum equal to S t , the number of strictly positive claim amounts is a binomial random variable with exponent t and parameter 1 -p 0 . For k ≥ 1, let S c k = X c 1 + . . . + X c k where the random variables X c i are i.i.d. and distributed as (X 1 | X 1 ≥ 1), i.e. P (X c i = n) = p n /(1 -p 0 ), n ≥ 1; put S c 0 = 0. By the law of total probability, we can then write that

P (S t = n) = p t 0 n k=0 t k 1 -p 0 p 0 k P (S c k = n), n ∈ IN, (2.19)
From (2.19), we thus see that the law of S t is of the form

P (S t = n) = p t 0 e n (t), (2.20)
where e 0 (t) = 1 and for n ≥ 1, e n (t) represents some polynomial of degree n in t satisfying e n (0) = 0. So far, we were only concerned with values of t ∈ IN. Formally, the right-hand side of (2.19), (2.20) might be evaluated for any value of t, positive integer or not. So, considering its natural extension on IR, we introduce the following quantities:

p t 0 e n (t
P (S t 1 +t 2 = n) = n i=0 P (S t 1 = n -i) P (S t 2 = i), n ∈ IN, (2.22)
and by inserting (2.20), this yields a formula of binomial type:

e n (t 1 + t 2 ) = n i=0
e n-i (t

= n) = n i=0 P (S t 1 = n -i) P (S t 2 = i), t 1 , t 2 ∈ IR. (2.24)
In the sequel, we will have recourse to the pseudo-probability mass functions of S t and S -t when t ∈ IN. These pseudo-distributions can be calculated numerically by recursion. For S t , it suffices to write that for t ≥ 0,

P (S t+1 = n) = n i=0 p n-i P (S t = i), n ∈ IN, (2.25)
with S 0 = 0. A similar procedure can be followed for S -t . Indeed, by (2.24) we have, for t ≥ 1,

P (S -t-1 = n) = n i=0 P (S -1 = n -i) P (S -t = i), n ∈ IN, (2.26)
so that it remains to find the pseudo-distribution of S -1 . For this, we recall that P (S 0 = n) = δ n,0 (the Kronecker delta), so that (2.24) with t 1 = -1 and t 2 = 1 yields the recursion

δ n,0 = n i=0 P (S -1 = n -i) p i , n ∈ IN. ( 2

.27)

Implications for T u . For clarity, we denote P (S t = n) = p * t n when negative values for t are allowed. We establish below that thanks to the previous pseudo-distributions, formula (2.16) can be generalized to the case where u ∈ IN.

Lemma 2.5 For any integers u ≥ 0, t ≥ 1 and u ≤ n ≤ u + t -1,

p * t n = n j=0 t + u -n t + u -j p * (t+u-j) n-j p * (j-u) j .
(2.28)

Proof. By (2.24) we can write that

P (S t = n) = n k=0 P (S t+u = n -k) P (S -u = k). (2.29)
Now, consider the event (S t+u = n -k) in (2.29). We first draw a straight line of slope 1 passing at point (t + u, n) (of equation y = -(t + u -n) + x); see Figure 4. Remember that t + u -n ≥ 1 by assumption. Thus, to reach the level n -k at time t + u, the process S i , 1 ≤ i ≤ t + u, necessarily crosses that straight line. Suppose that the first-crossing arises at some point (t + u -j, n -j) say, with k ≤ j ≤ n. Since this first-crossing time corresponding to τ t+u-n , we obtain, using (2.3) for P (τ t+u-n = n -j), that

P (S t+u = n -k) = n j=k t + u -n t + u -j p * (t+u-j) n-j p * j j-k . (2.30)
Inserting (2.30) in (2.29) then yields, after permutation of the two sums,

p * t n = n j=0 t + u -n t + u -j p * (t+u-j) n-j j k=0 p * j j-k p * (-u) k , (2.31)
which reduces to (2.28) since by (2.24), the sum ( j k=0 . . .) in (2.31) gives p * (j-u) j .

⋄

Note that p * (j-u) j in (2.28) is a standard probabilistic convolution when j ≥ u but this is not true for the first u values of j. By substituting (2.28) for p * t n in (2.12), we are able to deduce a new representation for the non-crossing probability of interest.

Lemma 2.6 For any integers u ≥ 0, t ≥ 1 and u ≤ n ≤ u + t -1,

P (S i < u + i for 1 ≤ i < t, and S t = n) = u j=0 t + u -n t + u -j p * (t+u-j) n-j p * (j-u) j . (2.32)
A remarkable property of formula (2.32) is that the sum is taken over u + 1 terms only, whatever the value of n. From that point of view, (2.32) is simpler than (2.12) where the sum contains n -u terms. Let us also point out that when u = 0, (2.32) directly reduces to the ballot formula (2.3).

For the non-ruin probability before t, inserting (2.32) in (2.17) provides the following alternative expression. A result of similar nature was obtained by Picard and Lefèvre [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] for a continuous-time version of the model. Proposition 2.7 (P.L.-type formula). P (T u = 1) = P (X 1 ≥ u + 1), and

P (T u ≥ t + 1) = u j=0 p * t j + p * (j-u) j u+t-1 n=u+1 t + u -n t + u -j p * (t+u-j) n-j , t ≥ 1.
(2.33) Formula (2.33) is given explicitly in terms of the pseudo-distributions of S -u , . . . , S t+u evaluated at different points between 0 and u + t-1. Comparing the efficacy of (2.18) and (2.33) is rather delicate (see, e.g., Rullière and Loisel [START_REF] Rullière | Another look at the Picard-Lefèvre formula for finite-time ruin probabilities[END_REF]). Some numerical experiments show that (2.33) becomes especially fast when u is small compared to t.

Infinite-time horizon

Let µ = E(X 1 ) be the expected amount of a claim. Since S t /t → µ as t → ∞, ruin will occur a.s. if µ ≥ 1. Assuming now that µ < 1, we want to obtain the probability of ultimate (non-)ruin. Specifically, from (2.18) and (2.33) we are going to establish the two formulas below.

Corollary 2.8 If µ < 1, then

P (T u < ∞) = (1 -µ) ∞ j=u+1 p * (j-u) j , (2.34)
P (T u = ∞) = (1 -µ) u j=0 p * (j-u) j . ( 2

.35)

Proof. To begin with, define by Z u+t-j a the expectation

Z u+t-j a = E 1 - S t+u-j t + u -j I (a≤S t+u-j ≤t+u-j) , (2.36)
for any naturals a ∈ [0, t + u -j] (I A is the indicator of the event A). Applying to (2.36) the dominated convergence theorem, we get that

Z t+u-j a → 1 -µ as t → ∞.
(2.37)

Now, consider the Seal-type formula (2.18). It contains a sum ( t+u+1 n=j . . .) that corresponds to Z t+u-j 0 . We then find that

P (T u < ∞) = 1 -lim t→∞ P (T u ≥ t + 1) = 1 -lim t→∞ P (S t ≤ u + t -1) + lim t→∞ u+t-1 j=u+1 p * (j-u) j Z t+u-j 0 ,
which reduces to formula (2.34) by virtue of (2.37) and since P (S t ≤ u + t -1) → 1.

For the P.L.-type formula (2.33), we recognize the inner sum ( u+t-1 n=u+1 . . .) as Z t+u-j u+1-j . By (2.37) and since P (S t ≤ j) → 0, we then obtain formula (2.35). ⋄ It can be checked that (2.34) and (2.35) correspond to formulas (40) and (42) obtained by Gerber [START_REF] Gerber | Mathematical fun with the compound binomial process[END_REF] for the usual version of the model with at most one claim per period. Let us point out that in practice, the sum (2.35) with u + 1 terms will be much easier to handle than the series (2.34).

Compound Poisson risk model

A continuous-time analogue of the previous model is the compound Poisson risk model. This model plays a prominent role in risk theory and a considerable literature is devoted to its study. The reader is referred to the books by Gerber [START_REF] Gerber | An Introduction to Mathematical Risk Theory[END_REF], Grandell [START_REF] Grandell | Aspects of Risk theory[END_REF], Panjer and Willmot [START_REF] Panjer | Insurance Risk Models[END_REF], Asmussen [START_REF] Asmussen | Ruin Probabilities[END_REF] and Kaas et al. [START_REF] Kaas | Modern Actuarial Risk Theory[END_REF]); see also, e.g., Gerber and Shiu [START_REF] Grimmett | Probability and Stochastic Processes[END_REF], Albrecher et al. [START_REF] Albrecher | On a gamma series expansion for the time-dependent probability of collective ruin[END_REF] and Cardoso and Waters [START_REF] Cardoso | Recursive calculation of finite time ruin probabilities under interest force[END_REF] for recent results on the ruin problem.

The insurance company is now assumed to be able to follow the evolution of its reserves in continuous-time. The occurrence of claims is described by a Poisson process N(t), t ≥ 0, with rate λ > 0. The successive claim amounts X i , i = 1, 2, . . ., are positive i.i.d. continuous random variables, with density function f (x) for x > 0. The initial reserves are of amount u ≥ 0, and the premium rate is constant and equal to c > 0. Thus, the reserves at time t > 0 are given by R t = u + ct -S t (as in (2.1)) but with S t = X 1 + . . .+ X N (t) and c > 0. Note that the law of S t has an atom at 0 with P (S t = 0) = e -λt , and is continuous on (0, ∞).

The ruin time is the first instant T u > 0 where R t becomes negative or null. In the sequel, we suppose that c = 1. This assumption, however, is not restrictive for the present model since it comes to change the monetary unit and claim amounts are here continuous. So, for t > 0, T u ≥ t means S y < u + y for 0 < y < t.

(3.1)
To obtain the exact distribution of T u , we will adapt the method of proof followed for the discrete-time model.

A standard approach

Basic case u = 0. The starting point is the following well-known ballot theorem in continuous-time (see, e.g., Takács [START_REF] Takács | Combinatorial Methods in the Theory of Stochastic Processes[END_REF]).

Lemma 3.1 (Ballot theorem, continuous case). Let S y , y ≥ 0, be a compound Poisson process with positive continuous jumps. Then, for any t > 0 and 0 < x < t, P [S y < y for 0 < y < t, and

S t ∈ (x, x + dx)] = t -x t f t (x)dx, (3.2)
where f t (.) denotes the density of S t on (0, ∞).

Let us recall that this result can be extended to processes with cyclically interchangeable increments. Various proofs of (3.2) can be found in the references mentioned earlier. We derive below formula (3.2) by arguing exactly as for (2.3).

Proof of (3.2). Denote by g t (x)dx the probability in the left-hand side of (3.2). Making a rotation of 180 • , g t (x) can be viewed as the density of the first-meeting time τ t-x of the process t -x + S y , y > 0, with the diagonal line of slope 1. In other words,

g t (x)dx = P [τ t-x ∈ (t, t + dx)]
= P [t -x + S y > y for 0 < y < t, and t -x

+ S t ∈ (t, t + dx)]. (3.3)
Considering the possible jumps during the time interval (0, t -x), we find that

g t (x) = x 0 f t-x (y) g x (x -y)dy. (3.4)
On another hand, since the increments of the process are independent and stationary, the following identity holds:

E(S t-x | S t = x) = (t -x) x t = x 0 y f t-x (y) f x (x -y) f t (x) dy. (3.5)
Let us define φ t (x) ≡ [(t -x)/t)] f t (x), for 0 < x < t; from (3.5) we see that these quantities satisfy the integral equation

φ t (x) = x 0 f t-x (y) φ x (x -y)dy. (3.6)
Comparing (3.4) and (3.6), we deduce that g t (x) = φ t (x), i.e. (3.2). ⋄

Turning to the ruin time T 0 , we can write that

P (T 0 ≥ t) = e -λt + t 0 P [S y < y for 0 < y ≤ t, and S t ∈ (x, x + dx)],
and thanks to (3.2), this yields formula (3.7) below, of structure similar to (2.10) (and going back to Cramér [START_REF] Cramér | On the Mathematical Theory of Risk[END_REF]).

Proposition 3.2 (Takács-type formula)

P (T 0 ≥ t) = e -λt + 1 t t 0 (t -x) f t (x)dx, t > 0. (3.7)
Passage to u ≥ 0. In the framework of the ballot problem, let us now suppose that the upper boundary is a straight line of slope 1 passing at point (0, u), u ≥ 0. Lemma 3.3 For any real numbers u ≥ 0, t > 0 and u < x < u + t, (1/dx) P [S y < u + y for 0 < y < t, and

S t ∈ (x, x + dx)] = f t (x) -e -λ(t+u-x) f x-u (x) - x u t + u -x t + u -z f t+u-z (x -z) f z-u (z)dz. (3.8)
Proof. After a rotation of 180 • , the probability in the left-hand side of (3.8), denoted by h t (x, u)dx, may be represented as h t (x, u)dx = P [τ t+u-x > t and S t ∈ (x, x + dx)].

(3.9)

Proceeding as for (2.13), we observe that to reach the level x at time t, the process t + u -x + S y , 0 < y ≤ t, will either cross or not the diagonal line, and if yes, at some level t + u -z where z ∈ (u, x]. Note that t + u -x is a possible crossing-level with the probability e -λ(t+u-x) . Thus, we get

f t (x) = h t (x, u) + e -λ(t+u-x) f x-u (x) + x u g t+u-z (x -z) f z-u (z)dz, (3.10)
where g t (x) is the density defined through (3.3). By the ballot formula (3.2),

g t+u-z (x -z) = t + u -x t + u -z f t+u-z (x -z), (3.11)
so that substituting (3.11) in (3.10) and using (3.9) yields formula (3.8). ⋄ Now, concerning the ruin time T u , we have, for t > 0,

P (T u ≥ t) = e -λt + u 0 f t (z)dz + u+t u
P [S y < u + y for 0 < y < t, and S t ∈ (x, x + dx)].

(3.12) Inserting (3.8) in (3.12) and permuting the two integrals, we then deduce the following result (originating from Prabhu [START_REF] Prabhu | On the ruin problem of collective risk theory[END_REF]). Proposition 3.4 (Seal-type formula)

P (T u ≥ t) = e -λt + u+t 0 f t (z)dz - u+t u e -λ(t+u-z) f z-u (z)dz - u+t z=u f z-u (z) u+t x=z t + u -x t + u -z f t+u-z (x -z)dx dz, t > 0. (3.13)
This formula is similar to (2.18) and shows how to obtain the probability of non-ruin before time t from the densities of S y where 0 < y < u + t.

An alternative approach

Pseudo-distributions for S t . First, let us consider the law of S t , t > 0. Apart from an atom at 0, it has a density on (0, ∞) given by

f t (x) = e -λt ∞ i=1 (λt) i i! f * i (x), x > 0, (3.14)
f * i (.) being the i-th convolution of the density of X 1 . Thus, similarly to (2.20), we can write f t (x) under the form

f t (x) = e -λt e x (t), (3.15)
where e x (t) represents the corresponding power series in t.

For t = 0, let e x (0) = 0. Now, let us extend the domain of t to IR in (3.15). It is directly checked that the series defining e x (t) converges absolutely. Therefore, as with (2.21), we may introduce a pseudo-density for S t by putting e -λt e x (t) ≡ ft (x) say, with t ∈ IR.

(3.16)

For t 1 , t 2 > 0, S t 1 +t 2 is distributed as S t 1 + S t 2 . In terms of the previous power series, this implies that e x (t 1 + t 2 ) = e x (t 1 ) + e x (t 2 ) +

x 0 e x-y (t 1 ) e y (t 2 )dy, x > 0.

(3.17)

Notice that (3.17) is an identity in t 1 (and t 2 ) that remains valid when t 1 < 0 (and t 2 < 0). Thus, the convolution property is still true for pseudo-densities, yielding

ft 1 +t 2 (x) = ft 1 (x) e -λt 2 + ft 2 (x) e -λt 1 + x 0 ft 1 (x -y) ft 2 (y)dy, t 1 , t 2 ∈ IR. (3.18)
These pseudo-densities can be rather easily calculated. Indeed, the number of claims at t > 0 being Poisson distributed, we have for f t (.), t > 0, the classical Panjer representation, namely

f t (x) = λt x 0 y x f (y) f t (x -y)dy, x > 0. (3.19)
Furthermore, it is clear that formula (3.19) remains true on IR, i.e. for ft (.). Now, (3.19) corresponds to a Volterra integral equation of the second kind, and such an equation can be solved numerically by appropriate discretization (see, e. g., Appendix D in Panjer and Willmot [START_REF] Panjer | Insurance Risk Models[END_REF]).

Implications for the law of T u . We will first derive a continuous variant of the identity (2.28).

Lemma 3.5 For any real numbers

u ≥ 0, t > 0 and u ≤ x ≤ u + t, f t (x) = e -λ(t+u-x) f x-u (x) + t + u -x t + u f t+u (x) + x 0 t + u -x t + u -y f t+u-y (x -y) fy-u (y)dy. (3.20)
Proof. We write f t = f t+u-u , so that applying (3.18) gives

f t (x) = f t+u (x) e λu + f-u (x) e -λ(t+u) + x 0 f t+u (x -y) f-u (y)dy. (3.21)
To begin with, let us consider the density f t+u (x -y) and a straight line of slope 1 passing at point (t + u, x). The compound Poisson process will reach the level x -y at time t + u necessarily after crossing that straight line, for the first time at some level x -z where z ∈ [y, x] (at time t + u -z). Notice that both levels 0 and x -y are possible first-crossing levels with positive probabilities. Using the density g t (x) defined by (3.3), we then get f t+u (x -y) = e -λ(t+u-x) f x (x -y) + g t+u-y (x -y) e -λy Note that (3.28) is constructed from the (pseudo-)densities fs (.), -u < s < t + u, evaluated at points between 0 and t + u. The same formula was obtained by De Vylder and Goovaerts [START_REF] De | Explicit finite-time and infinite-time ruin probabilities in the continuous case[END_REF], using a much less direct method based on the solution of some integrodifferential equation.

Infinite-time horizon

Let µ = E(X 1 ) and assume that λµ < 1. Then, ruin is not a.s. and the ultimate (non-)ruin probability can be expressed as follows. Proof. We observe that given any non-negative real a ∈ [0, t + u -y], the expectation

Z t+u-y a ≡ E 1 - S t+u-y t + u -y I (a≤S t+u-y ≤t+u-y) (3.31)
converges to 1 -λµ as t → ∞. Applying this result inside (3.13) and (3.28) then yields (3.29) and (3.30), respectively. ⋄ It can be checked that (3.29) is equivalent to the classical Pollaczeck-Khinchine (Beekman) formula. Formula (3.30), less standard, is given in De Vylder and Goovaerts [START_REF] De | Explicit finite-time and infinite-time ruin probabilities in the continuous case[END_REF]; see also Shiu [START_REF] Shiu | Calculation of the probability of eventual ruin by Beekman's convolution series[END_REF] for discrete claim amounts as discussed below.

Discrete claim amounts

Picard and Lefèvre [START_REF] Picard | The probability of ruin in finite time with discrete claim size distribution[END_REF] derived a similar formula (recalled in (3.37)) for a compound Poisson risk model where the claim amounts are positive integer-valued instead of continuous as before. The method of proof followed in that paper is primarily of algebraic nature. Let us sketch that the approach developed here, of more probabilistic essence, works well too.

Let {p t (n), n ∈ IN} be the probability mass function of S t at time t > 0, i.e. The following results can then be established. Suppose (for clarity) that u is a nonnegative integer, and denote by [t] the integer part of t. (Takács-type formula):

p t (n) = e -λt n i=1 (λt) i i! p * i n , (3.32
P (T 0 ≥ t) = 1 t [t] n=0 (t -n) p t (n), t > 0, (3.35)
(Seal-type formula):

P (T u ≥ t) = u+[t] j=0 p t (j) - u+[t] j=u+1 p j-u (j)   u+[t] n=j t + u -n t + u -j p t+u-j (n -j)   , t > 0, (3.36)
(P.L.-type formula):

P (T u ≥ t) = u j=0   p t (j) + pj-u (j) u+[t] n=u+1 t + u -n t + u -j p t+u-j (n -j)   , t > 0. (3.37)
These formulas are very similar to those derived for the compound binomial model (see (2.11), (2.18) and (2.33)).

A brief illustration

The Solvency II project provides a framework for a future prudential regulation of insurance in the European Union. In its Pillar 1, a main objective is to define the capital levels a company will have to allocate in order to cover all the risks inherent to the business. Two capital requirements are considered, the MCR (Minimum Capital Requirement) and the SCR (Solvency Capital Requirement) which has the key role. To determine the SCR, companies may follow standard approaches or apply internal models (the latter are more time-consuming but are expected to yield lower requirements).

A rather well accepted rule for the SCR consists in fixing the probability of solvency in any one year, most likely at the level of 99.5%. In other words, the SCR is such that the one-year 99.5% Value-at-Risk (V aR) is controlled. In recent discussions (Quantitative Impact Study 3), it has been suggested that the level of SCR could guarantee the MCR level to be available at anytime within 5 or 10 years, with a sufficiently high probability. Such a proposition would imply to deal with non-ruin probabilities in continuous-time over some given period (as investigated in the present paper).

Let us briefly illustrate that a finite-time non-ruin probability could indeed be used as an alternative to a Value-at-Risk. For that, we consider a compound Poisson risk model in which the Poisson parameter is λ = 1, claim amounts are exponentially distributed with mean µ = 1 and the premium rate is c = 1.1. First, we calculated the probabilities P (R t ≤ 0) of getting negative reserves at times t = 1, 5, 10 years, given initial reserves u = 0, 1, 2, 5, 10, 20, by using the Panjer algorithm; we also calculated the ruin probabilities ψ(u, t) over the period (0, t), for the same values of t and u, by applying the P.L.-type formula (3.37), as well as ψ(u) ≡ ψ(u, ∞) by using the well-known exact formula for the exponential case. The results are presented in Table 1. We observe that as expected, these probabilities can be very different, especially over a long horizon or for small initial reserves. Then, choosing several solvency thresholds 90, 95, 97.5, 99.5%, we proceeded in a similar way to determine the amounts u of initial reserves that enable to meet these thresholds, the criterion being either the Value-at-Risk V aR t or the non-ruin probability φ(u, t) with t = 1, 5, 10 years. The results, given in Table 2, show that the differences between required capitals can be important over a long horizon or for low solvency thresholds. Moreover, we note that, for example, an amount u = 6.03 needed to get a 1-year 99.5% Value-at-Risk allows also to get a 5-year non-ruin probability of ∼ 93%. In practice, however, controlling such a level of finite-time non-ruin probabilities is often found to be more reliable (see Loisel et al. [START_REF] Loisel | Robustness analysis, convergence of empirical finite-time ruin probabilities and estimation of risk solvency margin[END_REF]). So, in this framework, the choice of the 1-year 99.5% Value-at-Risk as the solvency criterion of reference could be questioned. A more detailed discussion on this problem is the object of a paper under preparation. u P (R 1 ≤ 0) ψ(u, 1) P (R 5 ≤ 0) ψ(u, 5) P (R

Figure 3 :

 3 Figure 3: First-crossing of {t -n + S i } through the diagonal line y = x.

Figure 4 :

 4 Figure 4: A trajectory {S i } with respect to the line of equation y = -(t u -n) + x.

g 0 f 0 e 0 g 0 f 0 e 0 g

 000000 t+u-z (x -z) f z (z -y)dz.(3.22) Now, in the right-hand side of(3.21), let us first examine the term f t+u (x) e λu . From(3.22) with y = 0, we see thatf t+u (x) e λu = e -λ(t-x f x (x) + e λu g t+u (x) + e λu x y g t+u-z (x -z) f z (z)dz. (3.23)Then, consider the integral in(3.21). Substituting(3.22) for f t+u (x -y), this integral becomesx t+u (x -y) f-u (y)dy = e -λ(t+u-x) x 0 f x (x -y) f-u (y)dy + x -λy g t+u-y (x -y) f-u (y)dy + s t+u-y (x -y) x 0 f z (z -y) f-u (z)dz dy.(3.24)By (3.17), the first integral in the right-hand side of (3.24) can be expressed asx 0 f x (x -y) f-u (y)dy = f x-u (x) -e λu f x (x) -e -λx f-u (x); (3.25)a similar expression holds for the third integral too. Making these two substitutions, we obtain for (3.24)x t+u (x -y) f-u (y)dy = e -λ(t+u-x) [f x-u (x) -e λu f x (x) -e -λx f-u (x)] + x -λy g t+u-y (x -y) f-u (y)dy + x t+u-y (x -y) [f y-u(y) -e λu f y (y) -e -λy f-u (y)]dy. (3.26) Finally, inserting (3.23) and (3.26) in (3.21), and using the ballot formula (3.3) leads to the result (3.20). ⋄ Lemma 3.5 allows us to obtain a new representation for the probability (3.8) of noncrossing before time t.

Lemma

Corollary 3 . 8

 38 If λµ < 1, then P (T u < ∞) = (1 -λµ) ∞ u f y-u (y) e λ(u-y) dy, (3.29) P (T u = ∞) = (1 -λµ) e λu + u 0 fy-u (y) e λ(u-y) dy . (3.30)

  IN, can be viewed as a pseudo-probability mass function for S t . The process S t , t ∈ IN, having independent stationary increments, the convolution property gives, for t 1 , t 2 ∈ IN,

) ≡ P (S t = n) say, with t ∈ IR.

(2.21)

Let us underline that the notation P (S t = n) in (2.21) has no probabilistic meaning except when t ∈ IN as in

(2.19)

. For example, when t = -1, P (S t = n) < 0 for all n ≥ 1. By analogy, however, P (S t = n), n ∈

  1 ) e i (t 2 ). (2.23) Now, we observe that the identity (2.23) still holds true for any t 1 , t 2 ∈ IR since both sides are polynomials of degree n in t 1 (and t 2 ) that are identical for an infinite number of values ot t 1 (and t 2 ). Therefore, (2.22) can be generalized to the pseudo-distributions above, giving P (S t 1 +t 2

  3.6 For any real numbers u ≥ 0, t > 0 and u ≤ x ≤ u + t, (1/dx) P [S y < u + y for 0 < y < t, and S t ∈ (x, x + dx)]

	=	t + u -x t + u	f t+u (x) +	0	u	t + u -x t + u -y	f t+u-y (x -y) fy-u (y)dy.	(3.27)
	Inserting (3.27) in (3.12), we then deduce the desired formula.
	Proposition 3.7 (P.L.-type formula)
	P (T u ≥ t) = e -λt +	0	u	f t (y)dy +	u	t+u	t + u -y t + u	f t+u (y)dy
	+	u y=0	fy-u (y)		u+t x=u	t + u -x t + u -y	f t+u-y (x -y)dx dy, t > 0.	(3.28)

)

  {p * i n , n ∈ IN} being the i-th convolution of the law of X 1 . As in (2.20), we write it as p t (n) = e -λt e n (t), (3.33) where e 0 (t) = 1 and e n (t), n ≥ 1, is a polynomial of degree n in t with e n (0) = 0. Then, for all t ∈ IR, we may define the pseudo-probability mass function {p t (n), n ∈ IN} of S t just by natural extension of (3.32), (3.33) to IR. In practice, these pseudo-distributions are calculated by the Panjer algorithm: for t ∈ IR,

	pt (0) = e -λt , and pt (n) =	λt n	n i=1	i p i pt (n -i), n ≥ 1.	(3.34)

Table 1 :

 1 10 ≤ 0) ψ(u, 10) ψ(u) .10 -4 3.1.10 -4 7.8.10 -6 9.2.10 -31.0.10 -7 0.032 0.366 20 8.6.10 -8 9.9.10 -8 6.0.10 -9 3.3.10 -5 6.0.10 -8 4.0.10 -4 0.148 Probabilities P (R t ≤ 0) of getting negative reserves at time t and ruin probabilities ψ(u, t) over the period (0, t), calculated for several values of t and u.Level V aR 1 φ(u, 1) V aR 5 φ(u, 5) V aR 10 φ(u, 10) φ(u)

	0	0.327	0.463		0.016	0.720	2.7.10 -4	0.785 0.909
	1	0.173	0.238		8.1.10 -3	0.512	1.2.10 -4	0.613 0.83
	2	0.089	0.120		3.8.10 -3	0.354	5.7.10 -5	0.470 0.758
	5	0.011	0.014		4.0.10 -4	0.103	5.3.10 -6	0.191 0.577
	10	2.599.5%	6.03	6.37 10.33	11.17	13.28	14.50 57.23
		97.5%	3.81	4.19	7.01	8.02	9.15	10.62 39.53
		95%	2.82	3.24	5.46	6.58	7.18	8.82 31.90
		90%	1.81	2.26	3.80	5.06	5.04	6.91 24.28

Table 2 :

 2 Initial reserves amounts u needed to satisfy different solvency thresholds, determined on the basis of the value-at-risk levels V aR t or the non-ruin probabilities φ(u, t), for several values of t.

Acknowledgements

We are grateful to the Referee for useful remarks and suggestions. This work was partly conducted while Cl. Lefèvre was visiting the Institut de Science Financière et d'Assurances, Université Claude Bernard Lyon 1. His sincere thanks go to all the members of the I.S.F.A. for their warm welcome and many interesting discussions.


L Gallimard

J P Pelle

Adaptive analysis in elastoplasticity based on an enhanced error in constitutive relation

Keywords: Finite element analysis, Error cause removal, Plasticity

Many industrial analyses require the resolution of complex nonlinear problems. For such calculations, error-controlled adaptive strategies must be used to improve the quality of the results. In this paper, adaptive strategies for nonlinear calculations in plasticity based on an enhanced error on the constitutive relation are presented. We focus on the adaptivity of the mesh and of the time discretization. Table X. Plastic work during the loading (SA solution)

Introduction

In elastoplasticity, the quality of the finite element solution depends on the parameters describing the discretizations used: the mesh elements' sizes, the order of the elements, the size of the time steps, but also the parameters of the iterative techniques used to solve the nonlinear problems (particularly the stopping criteria for the iterations). In elastoplasticity, low-degree elements (first-or second-degree) are generally used. Moreover, when the stopping criteria are sufficiently precise, the error introduced by the iterative techniques is negligible in comparison with the errors due to the mesh and the time discretization (Ladeve Áze et al., 1999). For these reasons, improving the quality of the finite element analysis requires, in most cases, simultaneous adaptations of the mesh and of the time discretization [START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF]. A classical approach is to perform an initial finite element analysis and to calculate a posteriori error estimators and error indicators on this first analysis. These errors are then used to determine new computation parameters in order to improve the quality of the calculation. The efficiency of the adaptive procedure depends greatly on the accuracy of the information given by the error estimator, both on the global level and on the local level. Three main approaches to the development of estimators can be found in the literature: the estimators introduced by Babus Ïka and Rheinboldt (1978), which use the equilibrium residuals to calculate the errors [START_REF] Babus Ïka | Computational error estimates and adaptive processes for some nonlinear structural problems[END_REF][START_REF] Johnson | Adaptative finite element methods in computational mechanics[END_REF][START_REF] Huerta | Error estimation and adaptative finite element analysis of softening solids[END_REF][START_REF] Rannacher | A a posteriori error estimation and mesh adaptation for finite element models in elasto-plasticity[END_REF]; the estimators introduced by Zienkiewicz [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical enineering analysis[END_REF], which consist of comparing the finite element solution with a smoothed solution [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical enineering analysis[END_REF][START_REF] Coupez | Adaptive solutions in industrial forming process simulation[END_REF]; the estimators introduced by Ladeve Áze [START_REF] Ladeve Áze | Me Âthode de calcul par encadrement des fre Âquences propres de structures e Âlastiques[END_REF], which are based on the concept of error in the constitutive relation [START_REF] Ladeve Áze | Accuracy of elastoplastic and dynamic analysis[END_REF][START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF]1997[START_REF] Gallimard | An enhanced error estimator on constitutive relation for plasticity problems[END_REF][START_REF] Ladeve Áze | A posteriori constitutive relation error estimators for nonlinear finite element analysis and adaptive control[END_REF]Ladeve Áze et al., 1999). The error measures in the constitutive relation satisfies the following property: if the error measure is equal to zero then the approximate solution of the problem is equal to the exact one. That means that this error measure takes into account all discretization errors over the whole time interval.

The purpose of the paper is, for finite element analysis in elastoplasticity, to compare the efficiency of adaptive techniques based on different versions of [START_REF] Drucker | On the postulate of stability of material in the mechanics of continua[END_REF] error in the constitutive relation. Three versions are compared:

(1) A version, called standard, introduced in Ladeve Áze et al. (1986) and developed in [START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF], which uses a rigorous quasi-explicit construction of the admissible fields. (2) A version, called ``p+k'', which uses an approximate numerical construction of the admissible fields. This version is an extension to elastoplasticity of techniques used in heat transfer by Strouboulis in [START_REF] Strouboulis | Recent experiences with error estimation and adaptivity, Part 1: Review of error estimators for scalar elliptic problems[END_REF] and in linear elasticity in [START_REF] Coorevits | Control of analyses with isoparametric elements in 2D and 3D[END_REF].

(3) An enhanced version developed in [START_REF] Gallimard | An enhanced error estimator on constitutive relation for plasticity problems[END_REF], where we introduce a minimization of the error on a set of admissible stress fields which leads to an enhanced admissible stress field which is nearer the exact stress field that the admissible field obtained by the standard and the ``p+k'' version. This version leads to better global and local effectivity indexes than the previous two.

In the first section, we briefly recall the formulation of the elastoplasticity problem which serves as reference. Then, the formulation of the associated discretized problem is given. Drucker's error on the constitutive relation is recalled, along with the principle of its implementation: construction of balanced force densities and construction of admissible stresses. The enhanced version and its different variants are presented and compared. The last section is devoted to the adaptation of the finite element analysis. To separate, in the global error, the contribution due to the spatial discretization from the contribution due to the time discretization, error indicators are proposed. The examples of adaptation presented confirm the effectiveness of the enhanced error estimator introduced in [START_REF] Gallimard | An enhanced error estimator on constitutive relation for plasticity problems[END_REF].

Problem to be solved Let us consider a structure , with boundary @. Over the interval 0; T the structure is subjected to:

. a prescribed displacement on a part @ 1 of the boundary: u d M ; t;

. a force density given by @ 2 @ À @ 1 : F d M ; t;

. a density of body forces in :

f d M ; t.
The problem which describes the evolution of the structure under the assumption of small perturbations is: Find uM ; t and M; t defined on  0; T such that:

. kinematic constraints and initial conditions

uM ; t P U 0;T V t P 0; T u j @ 1 u d V M P u j t0 0 V X 1 .
equilibrium equations

P S 0;T Vt P 0; T Y Vu à P U 0 Tr"u à d f d T u à d @ 2 F d T u à dS @ 2 .
constitutive relation

Vt P 0; T Y VM P j t A" u j Y t 3
In that formulation:

. U 0;T designates the space of the displacement fields and S 0;T the space of the stress fields. Generally these spaces are chosen in order to impose the condition that the energy of the fields be bounded on at each time t P 0; T. U 0 is the trial space for the displacement:

U 0 fu à such that u P U and u j @ 1 0g

where U is similar to U 0;T for the fields which are independent of time;

.

A is an operator which depends on the material and charaterizes the constitutive relation.

Finite element approximation

The discretization of the problem to be solved is classical. In space, a finite element model is used. For a given t, the finite element approximation u h of the displacement field is such that:

u h P U h
where U h is the discretized space chosen for U . The time interval 0; T is discretized into n subintervals

Á i t iÀ1 ; t i , Á f0 t 0 < t 1 < F F F < t n Tg
Classically, under the assumption that the displacement field is linear over each time increment, the approximate nonlinear problem to be solved at each time step t i is:

(2) Find u i h ; i h such that:

u i h P U h nd u i h u d t i on @ 1 4 Tr i h "u h à d f d t i T u h à d @ 2 F d t i T u h à dS 5 for each u h à P U h;0 i h À iÀ1 h A i "u i h À u iÀ1 h 6
where U h;0 U 0 U h and A i is the discretized operator which charaterizes the constitutive relation.

The non linear problem (4), ( 5) and ( 6) is solved using a Newton-Raphson algorithm. Thus, an error related to the precision of the algorithm is introduced. In general, the pair u i h ; i h given by the finite element solvers satisfies (4) and (6) exactly and ( 5) approximately. It is easy to perform an additional linear resolution to build a stress i h which satisfies (5) exactly. Under the assumption that the loading is linear over each time increment, the finite element solution is easily completed on 0; T. This solution will be called u h ; h . u h ; h satisfies (4) and ( 5) exactly and ( 6) approximately.

Error on the constitutive relation

The principle of the error on the constitutive relation consists of dividing the equations of the problem to be solved into two groups (Ladeve Áze and [START_REF] Ladeve Áze | Error estimate procedure in the finite element method and application[END_REF]:

. the admissibility conditions: the kinematic constraints and initial conditions (1) and the equilibrium equations (2);

. the constitutive relation (3). An approximate solution s ad u KA ; SA for the problem to be solved must satisfy the admissibility conditions (u KA satisfies (1) and SA satisfies (2)). The quality of this approximate solution is quantified by the non-verification of the constitutive relation (3).

Drucker's error on the constitutive relation For models which satisfy Drucker's inequality [START_REF] Drucker | On the postulate of stability of material in the mechanics of continua[END_REF] strictly, the quality of an admissible pair s ad u KA ; SA is evaluated by the error measure e (Ladeve [START_REF] Ladeve Áze | Accuracy of elastoplastic and dynamic analysis[END_REF][START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF]

À KA T " SA À " KA dt 8 and " KAjt " u KAjt KAjt A " KAj Y t SAjt A " SAj Y t V X 9
Remark: To calculate " SAjt , it is necessary to invert the constitutive relation. We also define the contribution of a time interval 0; t to the global error: e 0;t sup P0;t ; s ad d 10 and a relative error:

" e F D D 2 2 T 0 Â T KA " KA T SA " SA Ã d dt 11
Construction of an admissible solution The difficulty is to derive an admissible solution s ad u KA ; SA from the finite element solution u h ; h . For the displacement, we simply choose u KA u h . For the stress, we follow the method described in Ladeve Áze et al. (1986) and [START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF], which uses the construction techniques developed in elasticity [START_REF] Ladeve Áze | Error estimation and mesh optimization for classical finite elements[END_REF]. We briefly recall the main ideas of this construction.

In a first step, we derive from the finite element solution i h M at t i a stress field i SA M which satisfies the equilibrium equation at t i exactly (2):

i SA P S Vu à P U 0 Tr i SA "u à d f d t i T u à d @ 2 f d t i T u à dS @ 12
where S is similar to S 0;T for the fields which are independent of time.

In a second step, we interpolate the balanced fields linearly over each time step. Under the assumption that the loading is linear over each time step (an assumption which, in practice, is not very restrictive), the constructed field is admissible at each instant.

The construction of i SA is performed in two substeps. Substep 1. We build on the edges of each element E of the mesh the traction forces F i h which represent the stress vectors on the edges of the element.

SAjE n E F i h on @E 13

(the function E depends on the element E; it can assume two constant values: +1 or ±1. Moreover, for two adjacent elements, we prescribe E H E 0 on the common edge À EE H .) Naturally, on the edges belonging to @ 2 , E F i h F d t i is prescribed. Moreover, the tractions F i h are constructed such that each element of the mesh is in equilibrium:

for each rigid body displacement u S E f d t i T u S dE @E E F i h T u S dÀ 0: 14
The key to the method is the prolongation condition which links i SA to the finite element solution:

For each element E in the mesh and for each node

j E i h À i SA T " j dE 0 15
where j is the shape function associated with node j.

The main lines of the construction of F i h are given in the Appendix. Substep 2. i SA is built on each element E from the element tractions F i h as the solution to the local problem:

div i SA f d t i 0 in E i SA n E F i h on @E @ 16
Two techniques may be used. The first one, described in Ladeve Áze et al. (1991), consists of determining a simple polynomial solution. More precisely, each triangular element is subdivided into three subtriangles built on the inertia center and the three vertices. We assume that on each subtriangle i SA is linear for linear elements and quadratic for quadratic elements. This construction is quasi-explicit and i SA is a strictly admissible stress field. This technique leads to the standard version of the estimator of the error on the constitutive relation. This error estimator will be designated by e std (the corresponding relative error is designated by " std ).

A second technique consists of constructing on E the stress field i SA as a solution to the following problem:

pind i SA tht stisfies IU suh tht i SA minimizes E TrK À1 dE & 17
where K is a definite symmetric positive operator. The dualisation of this problem leads to:

pind v i E displement field defined on E for eh v à defined on E X E TrK "v i E "v à dE E f t i T d v à dE @E F i h T v à dS i SAjE K "v i E 18 (14) ( 15
)
This problem is solved numerically on E by the finite element method. We use only one finite element of degree p k (where p is the degree of the initial finite element analysis and k a positive integer). This technique, which yields an approximate admissible field, was introduced in [START_REF] Strouboulis | Recent experiences with error estimation and adaptivity, Part 1: Review of error estimators for scalar elliptic problems[END_REF] for thermal problems. The corresponding version of the error is called p k and the estimator is designated by e pk . (The corresponding relative error is designated by " pk .) Remark: a priori this technique works with all definite symmetric positive operator K. In practical computations we use operators linked with the F.E. analysis, for instance the Hooke operator, the tangent rigidity at t i or at t iÀ1 . In practice the numerical examples show that the obtained error measure is little sensitive to this choice.

Enhanced Drucker constitutive error estimator Principles

An enhanced version of the Drucker constitutive error estimator, derived from a previous work in elasticity [START_REF] Ladeve Áze | New advances on a posteriori error on constitutive relation in f.e. analysis[END_REF], was presented in [START_REF] Gallimard | An enhanced error estimator on constitutive relation for plasticity problems[END_REF]. The aim is to optimize the construction of the admissible field SA in order to improve the quality of the associated error estimator.

The principle of the method may be described as follows:

.

The prolongation condition ( 15) is weakened by applying it only to the non-vertex nodes:

VE P T h nd for eh non-vertex node j E i h À i SA T " j dE 0 19
As a result, on each element edge À the traction F can be split into two parts:

F j À R j À H j À 20 H j À is a traction with zero resultant force and moment which is completely determined by the weak prolongation condition (19). R j À is a traction which depends on two arbitrary vectors R 1jÀ and R 2jÀ which need to satisfy only the global equilibrium (14).

. For each set of R j À which satisfies the equilibrium conditions, a stress field i SA R that belongs to S can be constructed using the techniques defined in the preceeding paragraph. R is the set of vectors R jjÀ j 1; 2. If we assume that enh;iÀ1 SA has already been calculated, R opt is determined by minimizing the error contribution on each time step t iÀ1 ; t i approximately:

J R t i t iÀ1 SA R À KA T " SA R À " KA dtd 21 where SA R t À t iÀ1 t i À t iÀ1 i SA R t i À t t i À t iÀ1 enh;iÀ1

SA

The enhanced balanced stress at t i is designated by enh;i SA i SA R. The calculated stress field on 0; T is designated by enh SA . The details of this technique can be found in [START_REF] Gallimard | An enhanced error estimator on constitutive relation for plasticity problems[END_REF], together with examples which show that a very good quality error estimator is obtained. In the following sections, this enhanced error estimator will be designated by e enh (the corresponding relative error will be designated by " enh ).

Variants to the enhanced estimator

In practice, the optimization is not performed over the whole set of element edges. It is sufficient to perform the minimization on the element edges belonging to ``sensitive'' zones, e.g. high stress gradient, transition zone between plastic and elastic regions, deformed elements, F F F In the other regions, we use the standard densities calculated by the method presented in the Appendix. The choice of the regions where the optimization is performed obviously plays an important role in the method.

Two other parameters may a priori have an influence on the quality of the stress field SA obtained:

(1) The number of integration points n Á used on Á i t iÀ1 ; t i to evaluate the function ( 21).

(2) The criterion chosen to stop the conjugate gradient algorithm used to solve the optimization problem.

The purpose of the following sections is to study the influence of these parameters on the quality of the associated error estimator.

Definition of the numerical tests

The numerical tests presented were performed on the structure shown in Figure 1. This structure is subjected to a prescribed displacement: ũd u 0 :t:x (where u 0 is a scalar quantity). The material parameters are:

. Young's modulus E 200:E3 Mpa, .

Poisson's ratio 0:3, .

Initial yield stress 0 240 Mpa, .

Isotropic hardening law R 500p 1=2 .

Two calculations were performed: the first on a quasi-regular mesh (Figure 2 (2,028 DoF) with six time increments), the second on a refined mesh (Figure 3 (9,672 DoF) with 16 time increments). The meshes were made of 6-node triangular elements. The Von Mises stresses at the end of the loading are shown on Figure 4.

The values of the error estimators e std and e pk (as well as the relative values " std and " pk ) for Calculations 1 and 2 are summarized in Table I.

Stopping criterion and influence of n Á

In the examples studied, the function J R decreases rapidly during the first iterations of the conjugate gradient technique. Since our objective is to make J R decrease and not to determine the best possible value of R, it is not necessary to iterate until convergence. The stopping criterion of the conjugate gradient technique is based on the mean variation of J R during two consecutive iterations.

The study of the influence of the stopping criterion was performed on Calculation 1 (Table II) (The optimization was performed on the whole mesh with n Á 8). A stopping criterion of 5 per cent seems sufficient to obtain a great improvement: the error decreased from e pk 1:09 to e enh 0:51. The influence of n Á for the calculation is shown in Table III. (The optimization was performed for the whole mesh with a stopping criterion equal to 1.10-4.) We observe that this parameter has very little influence on the error estimator.

Criterion for choosing the patches

Let us define the error density d E on an element E of the mesh T h :

d E e E mesE p 22
where mesE is the area of element E. Let d max be the maximum density:

d max sup EPT h d E 23
The simplest criterion to determine the ``rough'' zones and choose the edges À where R À must be optimized is to calculate the error e pk and to determine the elements where the error density is greater than a fraction of d max . The selected elements with common edges are gathered together into element patches. It is useful to add one or two layers of elements around each patch in order to increase the chance of minimizing the error contribution on elements with a high error density. Then, the patches with common interior edges are gathered together, yielding n patches ! k . The standard tractions are prescribed on the edge of each patch ! k ; then J j ! k is minimized. The default values considered are 0 (all elements are included in the function to be minimized), the stopping criterion of the minimization algorithm is 10 À2 and n is taken equal to 1.

Influence of the choice of patches

The results obtained for Calculation 1 with the error estimator e enh are shown in Table IV. We observe that, on this coarse mesh with a large plastic region, the improvement due to the error estimator e enh is significant, even for a coefficient 0:5 (e pk 1:09 and e enh 0:59). The patches chosen are on Figure 5).

On Calculation 2, which was performed on a refined mesh, we observe that for a coefficient 0:5 we again obtain a significant improvement in the error estimation (Table V). The patches used for 0:5 are shown on Figure 6.

Adaptive control technique

The adaptive technique proposed relies on an initial coarse calculation. From this first calculation, we determine the parameters (mesh, time discretization) of a second calculation in order to obtain a given precision " 0 for as little cost as possible. The error measure " is global in space and time and takes into account the mesh quality, the quality of the time discretization and the precision of the stopping criterion of the iterative technique used. This global information is not sufficient to determine the calculation parameters to obtain the required precision. To determine these parameters, we use tools which we call error indicators which allow us to evaluate the errors due to the mesh and the errors due to the time discretization.

Error indicators

We use the time error indicator developed in [START_REF] Gallimard | Error estimation and adaptivity in elastoplasticity[END_REF], which is common to the three global error estimators. The calculation of the space error indicator depends on the construction technique of the admissible stress field. We use the definition of the space error indicator given in Ladeve Áze and Moe Ès (1999). Basically, this indicator is constructed by choosing as the reference problem the problem which is discrete in time and continuous in space. The error occurring between this new reference problem and the finite element solution is due only to the space discretization.

An admissible solution for this problem is the set of pairs s i ad u i KA ; i SA for i P f1; F F F ; ng.

Let us define: where Áx i1 x i1 À x i 25 and

space s i ad i 1 SA À i1 KA T Á" i1 SA À Á" i1 KA 24
i1 À i A i " i1 À " i 26
A i is the discretized operator of the constitutive relation obtained by an implicit Euler integration scheme.

A new space error indicator may be defined with the enhanced admissible pairs s enh;i ad u i KA ; enh;i SA for i P f1; F F F ; ng:

I enh space 2 sup kPf0;FFF;nÀ1g k i0 space s enh;i1 ad d 27
It is possible, like for the global error, to define a relative space error indicator: i enh space . The space error indicators i enh space and i enh space are calculated in the same way.

Relation between ", i space and i time Let us consider the quantity i tot defined by:

i 2 tot i 2 space i 2 time 28
In Ladeve Áze and Moe Ès (1998), it is shown numerically that:

i std tot i std space 2 i 2 time q 9 " std : 29

The same property is observed for the ``enhanced'' and ``p+k'' quantities. On Table VI, we increase the quality of the mesh for a given time discretization. On Table VII, we increase the number of time steps for a given mesh.

Adaptation of the calculation parameters

In the following section, since the adaptive technique proposed is the same for the three error estimators, we will omit the exponents std, p k and enh. For a prescribed global error " 0 , the following procedure may be defined:

. An initial analysis over the time interval 0; T is performed on a coarse mesh T h for a time discretization Á. Then, the global error ", the time indicator i time and the space indicator i space are calculated.

.

A new mesh T Ã h and a new time discretization Á Ã are determined in order to minimize a cost function f N Ã ; M Ã M Ã :N Ã with the restriction:

i à space 2 i à time 2 " 2 0
where N Ã is the number of elements in T Ã h and M Ã is the number of time steps in Á Ã . This minimization problem is solved by using techniques introduced in Ladeve Áze et al. (1991) and later in [START_REF] Coorevits | Mesh optimization for problems with steep gradients[END_REF] and Ladeve Áze and Moe Ès (1999).

Remarks:

. For the 6-node triangular elements that we use in the examples, the value of is set to 1 and the value of is set to 2.

. If the prescribed global error " 0 is too small compared to the global error " obtained in the initial calculation, it is, in practice, difficult to determine the parameters of the optimized calculation directly. One or more intermediate calculations must be performed with target errors equal to " (where, for instance, 0:25).

Examples

To illustrate this procedure, we reexamine the example of the structure subjected to a prescribed horizontal displacement 1. The initial analysis (Figure 7) is performed with two time steps (this calculation is referred to as Calculation 0).

The three error estimators presented in the previous section are used in succession with the adaptive procedure. The parameters chosen for the calculation of " enh are: n Á 1, stopping criterion = 5:10 À 2 and 0:5.

The prescribed global error " 0 is 5 per cent.

In Table VIII, we show the error, the number of nodes and the number of time steps obtained for the optimized calculation. Figures 8-10 show the resulting discretizations.

A zoom is performed on Zone A to study the mesh more precisely. The error measure " std is very large in the presence of a high stress gradient and that leads to an overrefinement of Zone A (Figure 11). Moreover, the high error gradient in this zone leads to the creation by the mesh generator of deformed Optimized calculation e std ± mesh zone A elements. For the error estimator " pk , the overrefinement is, again, important (Figure 12). For the error " enh (Figure 13), the minimization performed in the zones where the error is overestimated enables this phenomenon to be corrected, resulting in a more balanced mesh. Thus, the size variations of the elements are more easily dealt with by the mesh generator.

The evolutions of some quantities of interest were also studied during the evolution of the adaptive procedure. Calculation 1 corresponds to the discretization obtained after one iteration of the adaptive technique. Calculation 2 corresponds to the discretization obtained when the global error reaches approximately 10 per cent.

In Table IX, we represented the evolution of the plastic work calculated from the finite element solution. The three calculations converge toward the same result.

In Table X we represented the evolution of the plastic work calculated from the statically admissible solution. The three calculations converge toward an identical result. But we observe that for calculations 0 and 1, which are very coarse, the admissible solution is not good and the improvement due to the optimization technique is very significant.

In Table XI, we show the evolution of the size of the plastic zone (calculated from the finite element solution). Once again, the three calculations converge toward similar values and the values calculated during the initial analysis are very inaccurate.

Conclusion

In this paper, we have proposed a study of the enhanced error " enh on the constitutive relation in plasticity and shown that the calculation of this error estimator in order to estimate the global quality of the calculation may be restricted to the zones which create difficulties. We have shown that the use of this error estimator and associated error indicators in an adaptive procedure allows us to end up with less refined and more regular meshes while preserving the precision and the quality of the calculation.

Appendix. Standard recovery method

The construction of the balanced stress field is performed at each time step. For the sake of simplicity, we will omit the exponent std; i.

Prolongation condition

The admissible stress SA is linked to the finite element stress field by the ``strong'' prolongation condition defined by the relation ( 15). This can be written using the equilibrium of SA :

@E SA :n T j dÀ E h T " j À E f T d j dE Q E j 30
This relation introduces the element tractions on the element boundaries. The element traction on element E is designated by E F h .

. F h is defined on the element edges.

. E is a scalar function defined on E such that E AE1 and E H E 0 on the common edge À EE H of two adjacent elements E and E H .

With the expression:

E F h SA :n on @E 31 Equation ( 30) becomes:

@E E F T h j dÀ E h T " j À E f T d j dE Q E j 32
Remark: On an edge À P @ 2 we impose E F h F d In Ladeve [START_REF] Ladeve Áze | Me Âthode de calcul par encadrement des fre Âquences propres de structures e Âlastiques[END_REF], it is shown that if E F h satisfies the prescribed loading on @ 2 and Equation (32), then the element tractions F h are balanced.

Standard construction of F h

The complete construction can be found in Ladeve Áze et al. (1991). We will consider only the case where the nodes are internal to .

Non-vertex node: j P À EE H common edge to E and E H . The prolongation condition (32) leads to:

À EE H E F T h j dÀ Q E j 33 À EE H H E F T h j dÀ Q H E j 34
Writing the finite element equilibrium leads to:

Q E j Q H E j 0 35
Therefore,

À EE H F T h j dÀ H E Q H E j E Q E j 36
Vertex node j: We write (cf Figure A1):

b k À k E k F T h j dÀ 37 Figure A1.

Vertex nodes

Equation (32) for the elements connected to node j leads to the following linear system:

b 1 À b 2 Q E 1 j b 2 À b 3 Q E 2 j b 3 À b 4 Q E3 j b n À b 1 Q En j V b b b b b b b b X 38
This system has an infinity of solutions if kn k1 Q E k j 0, a property which occurs if h satisfies the finite element equilibrium. In practice, b i are chosen in order to minimize: hjE k :n k À hjE kÀ1 :n kÀ1 ! j dÀ. The element tractions are then recovered by assuming that F h is a linear combination of the shape functions on the element edges:

F i h A i A i
jP " I A i j j 40

Figure 1 .

 1 Figure 1. Problem to be solved

  Figure 2. Bracket ± Calculation 1

  Figure 3. Bracket ± Calculation 2

  Figure 4. Bracket ± Calculation 1 ± Von Mises stress ± Initial yield stress 240Mpa

  Figure 7. Calculation 0

  Figure 9. e pk ± optimized calculation

  Figure 12.Optimized calculation e pk ± mesh zone A

f b 1

 1 ; b 2 ; :::;

  :

	e sup i 1=2 EC h t; s ad d where 0 t 18,7 t; s ad SA
	tP0;T

Table IV .

 IV Downloaded by Université Paris Ouest Nanterre La Défense At 05:43 08 January 2018 (PT)

	Coef.	0.5	0.4	0.3	0.2	0.1	0.0
	e enh e pk À e enh =e pk e std À e enh =e std	0.59 46% 55%	0.55 49% 58%	0.48 56% 64%	0.48 56% 64%	0.44 60% 67%	0.44 60% 67%	Influence of ± Calculation 1
								Figure 5.
								Calculation 1

Table V .

 V

		Coef.	0.5	0.4	0.3	0.2	0.1	0.05
	Influence of on Calculation 2	e enh e pk À e enh =e pk e std À e enh =e std	0.267 57.2% 58.4%	0.255 59.1% 60.3%	0.192 69.2% 70.1%	0.167 73.2% 74.0%	0.160 74.3% 75.1%	0.158 74.7% 75.4%
	Figure 6.

Table VIII .

 VIII Optimized calculation

		e std	e pk	e enh
	Error (%)	4.89	5.01	4.70
	NDDL	7,714	6,356	5,206
	Number of steps	7	7	7
				Figure 8. e std ± optimized
				calculation

Table IX .

 IX Plastic work during the loading (KA solution)

	Calculation	e std	e pk	e enh
	0	2.66	2.66	2.66
	1	2.51	2.49	2.51
	2	2.34	2.34	2.34
	3	2.33	2.33	2.33


Claude Tadonki
email: claude.tadonki@mines-paristech.fr

Setting up an Operational Academic HPC Center

Hardware support for high-performance computing (HPC) has so far been subject to significant advances. The pervasiveness of related systems, mainly made up with parallel computing units, makes it crucial to spread and vivify effective HPC curricula. Besides didactic considerations, it appears very important to implement HPC hardware infrastructures that will serves for practices, and also for scientific and industrial requests. The latter ensures a valuable connection with surrounding cutting-edge research activities in other topics (life sciences, physics, data mining, applied mathematics, finance, quantitative economy, engineering sciences, to name a few), and also with industrial entities and services providers from their requests related to HPC means and expertise. This aspect is very important as it makes an HPC Center becoming a social actor, while bringing real-life scenarios into the academic context. The current paper describes the major steps and objectives for a consistent HPC curriculum, with specific analyses of particular contexts; suggests how to technically set up operational HPC infrastructures; and discusses the connection with end-users, all these in both effective and prospective standpoints.

I. INTRODUCTION

With the advent and pervasiveness of multicore processors, parallel processing is continuously getting out of mystery. However, even if it is true that most of end-users of common applications are now somehow aware of their parallel execution, and that ordinary programmers are ready or more willing to admit the necessity of parallel processing considerations, delving into parallel computing topics [START_REF] Grama | Introduction to Parallel Computing (2nd Edition)[END_REF], [START_REF] Leighton | Introduction to Parallel Algorithm and Architectures: Arrays, Trees, and Hypercubes[END_REF], [START_REF] Lewis | Introduction to Parallel Computing[END_REF], [START_REF] Roosta | Parallel Processing and Parallel Algorithms: Theory and Computation[END_REF], [START_REF] Gebali | Algorithms and Parallel Computing[END_REF], [START_REF] Jaja | Introduction to Parallel Algorithms[END_REF], [START_REF] Leighton | Introduction to Parallel Algorithms and Architectures: Arrays, Trees, Hypercubes[END_REF] is still not an instinctive nor intuitive step. People need to be convinced about that choice, and the balance between the underlying effort and the corresponding reward should appear favorable. Thus the need of incentive HPC initiatives.

Significant efforts have been made so far to make parallel computing implementation as easy as possible, if not seamless from the programming viewpoint [START_REF] Karniadakis | Parallel Scientific Computing in C++ and MPI: A Seamless Approach to Parallel Algorithms and their Implementation[END_REF], [START_REF] Kirk | Programming Massively Parallel Processors: A Hands-on Approach[END_REF], [START_REF] Herlihy | The Art of Multiprocessor Programming[END_REF], [START_REF] Pacheco | Parallel Programming with MPI[END_REF], [START_REF] Quinn | Parallel Programming in C with MPI and OpenMP[END_REF]. However, the topic is still non trivial [START_REF] Correa | Models for Parallel and Distributed Computation: Theory, Algorithmic Techniques and Applications[END_REF], [START_REF] Hromkovic | Communication Complexity and Parallel Computing[END_REF]. To date, there is no effective compiler who can automatically derive a parallel version of a given ordinary sequential code without explicit user annotations or monitoring. Indeed, parallel computing adds another abstraction to basic algorithm design that has nothing to do with the underlying problem being solved. Thus, the topic might be perceived as a plain complication, otherwise technically speculative. Thus, we need to grab people attention and make them accept to take the way of parallelism. This is an education goal, prior to targeting an increasing number of (potential) candidates for parallel computing curriculum.

There is a increasing need of programmers with parallel programming skills. Thus, it is now common to see job ads seeking people qualified or specialists in parallel computing. This has significantly motivated parallel computing curricula, which, from our point of view, should not stand totally disconnected from other areas. Indeed, it is always better to understand what we have to parallelize. This will result in a more skillful parallel scheme and yield a consistent selfcontained report. Thus, it is preferable to have good basic skills before moving to parallelism. This through has to be kept in mind while designing an HPC curriculum.

From a technical point of view, especially with shared memory parallelism [START_REF] Chapman | Using OpenMP: Portable Shared Memory Parallel Programming[END_REF], there is a heavy hardware concurrency, which should be managed or taken into account carefully at design or evaluation time [START_REF] Rauber | Parallel Programming for Multicore and Cluster Systems[END_REF]. Thus, it is important to make sure that the basis of computer architecture [START_REF] Mano | Computer System Architecture[END_REF], [START_REF] Patterson | Computer Architecture: A Quantitative Approach Fifth Edition[END_REF] are well understood, maybe making it a prerequisite. Another reason for this is the impact that a parallel scheme might have on numerical aspects. Also for proof of correctness, hardware mechanisms of parallel systems should be included into the equation. Unfortunately, it is very common to encounter computer architecture unaware parallelizations. This might be due to a lack of related skills or just a plain carelessness. The most noticeable example is the important proportion of programmers who are completely unaware of vector computing, also referred as SIMD computing, whose related hardware elements and mechanisms are implemented in almost all modern processors. Thus, providing a computer architecture course at the earliest stage is certainly a good idea.

We see that there are numerous facts and thoughts that emerge when it comes to parallel computing from a didactic viewpoint. Thus, we need to efficiently organize and manage all related teaching/training/expertise activities [START_REF] Tadonki | Basic parallel and distributed computing curriculum[END_REF]. Beside invariant common needs, there are specific contextual issues that need to be addressed conscientiously through well targeted applications. Within an academic context, the best way to handle all these initiatives is to set up an HPC Center that will implement the main didactic actions and promote the use of parallel computing as previously stated. The purpose of this paper is to provide a consistent panoramic view of an academic HPC Center including the key topics with their interaction, and to describe its main associated activities. The rest of the paper is organized as follows. The next section describes key didactic factors related to HPC, followed by an overview of how to design and implemented an HPC Center. Section IV addresses the question of building a local HPC cluster. Section V concludes the paper.

II. DIDACTIC KEY FACTORS ABOUT HPC

A. Justify why considering parallel computing

First of all, we have to sell the need for HPC [START_REF] Buyya | High Performance Cluster Computing: Architectures and Systems[END_REF], [START_REF] Severance | High Performance Computing[END_REF], [START_REF]High Performance Computing in Science and Engineering[END_REF] to our (potential) audience. This preliminary communication is crucial, because people need to be convinced that the related investment is worth considering. The requested effort is not only technical, but also conceptual, as it implies a change of mind about computing devices, programs, algorithms and programming. This should be a clear winning step in order to capture a sincere and keen attention of the listeners, their enthusiasm will determine the success of the course and associated exchanges. We should always keep in mind that ordinary people expect the increasing CPU frequency to fit their needs for speed. The way to achieve the aforementioned marketing goal is to raise a number of (intentionally pro domo) arguments. We state some of them:

• There are numerous important real-life applications where processing speed is essential, but which cannot be fulfilled at the expected level by a single CPU [START_REF] Greenlaw | Limits to Parallel Computation: P-Completeness Theory[END_REF], [START_REF] Dongarra | Introduction to the HPC Challenge Benchmark Suite[END_REF], [START_REF] Tadonki | High Performance Computing as Combination of Machines and Methods and Programming[END_REF]. Nice and illustrative examples should be described here, some of them following an anecdotal way if possible. Complicated cases should be avoided, even if (and this can be at most mentioned) scientific computing is the major "customer" of parallel computing [START_REF] Hager | Introduction to High Performance Computing for Scientists and Engineers[END_REF], [?], [START_REF] Karniadakis | Parallel Scientific Computing in C++ and MPI: A Seamless Approach to Parallel Algorithms and their Implementation[END_REF]. We suggest for instance: urgent query in a large distributed database; real-time high definition image processing; realistic high precision simulation; new generation video gaming. With a scientific audience, we use to display figure Fig. 1 with difficult network graphs problems as examples. • The so-called Big Data is an active topic for good reasons. Indeed, the volume of data to be processed in number of real-life applications is more and more huge, and the response is expected to remain almost instantaneous whatever the scenario. The case of social networks is a good example here. Data analytics applications can be also mentioned. The major actors of the two previous topics are developing a significant research activity in parallel computing, with explicit laboratories implemented in various countries. • Almost all modern processors integrate several CPU cores. This fact by itself can also be considered as a plain argument for delving into parallel computing concepts. Even end-users of applications that run in parallel could get a better technical picture from being aware of that. Personally, we use to introduce the current fact by considering the OS multitasking, which allows to concurrently run several applications even with a single CPU core. Then explain that this virtual concomitance becomes effective with parallel units. • Energy is a major concern in most of modern activities including those related to the use of computers. It should be recall here that energy is nearly proportional to running time. Thus, reducing the energy consumed by a given software application can be achieved by accelerating its execution. In addition, a processor with clock frequency f consumes more energy than a dual-core with frequency f /2 per core. The slogan here is "parallel computing is energy efficient". • Many natural activities are intrinsically parallel. Thus, natural parallelism is intuitively expected to be so implemented on computers. Illustrative examples can be taken from simulation and artificial intelligence.

B. Keep HPC useful for others

This is very important in order to make the HPC Center being perceived as a useful entity and thereby receive rewarding considerations. In addition, keep addressing external requests related to industrial activities, scientific research, societal issues,services, and trainings, will certainly confirm the importance of the Center and vivify its activities. In practice, it is common to see financial sponsors firmly recommend such connection and also use it a crucial metric for their evaluation and decision. The same exhortation might come from institutional hierarchy or politics. Note that, because of the strong need for HPC availability, some companies implement their own HPC section. Academics could claim permanent investigations, less restricted research focus, and high potential for innovation.

C. Make sure to address local needs

From our personal teaching and collaboration experience, in South America for instance, it clearly appeared that addressing local needs and issues is highly appreciated and valuated. This is normal because otherwise, especially in case of too specific concerns, there is no chance that an appropriate solution will come from elsewhere. In addition, there is a strong desire to have direct hand to the implemented solutions for many reasons including cost, suitability, independence, and technical maturity. Thus, an HPC Center should prioritize local issues in their applications portfolio, of course with a certain of level interest in external and prospective scenarios.

D. Be keen in academic collaborations

Standing and acting in an isolated manner is a noticeable factor of stagnation and sort of obscurantism, which is the opposite of what is expected from an academic entity. HPC is moving too fast and changes occur very frequently. Hardware advances and near future projections are regularly announced. New HPC methods, achievements, and challenges are actively made available. Thus, keep collaborating with other actors is essential to survive and to receive consideration. The resulting synergy is another vital benefit to take into account.

III. IMPLEMENTING AN ACADEMIC HPC CENTER

This section discusses the main components and objectives of an academic HPC Center.

A. Global organization of the Center

The Center itself is a virtual entity of people with relevant skills. Figure Fig. 2 provides a view of the main units of the Center. Beside this logical clustering, it is important to have the units interacting efficiently. External requests as previously explained are likely to reach the head or research staff. A management policy should (at least roughly) indicate how process such requests.

B. A consistent HPC curriculum

We present and discuss the main basis for a consistent HPC curriculum [START_REF] Tadonki | Basic parallel and distributed computing curriculum[END_REF]. Still at the level of an overview, we provide a connected panorama of the key learning units related to parallel computing. The teaching part could be complemented with external invited speakers, thus the importance of collaborations, which will also boost the outcome of the research activities.

1) Topics selection: The curriculum should be organized so as to allow coherent inner steps and motivated choices. Figure Fig. 4 shows the main paradigms of parallel programming. These paradigms are independent, so the associated classes can be provided without any strong interconnection. However, it should be made possible to follow all of them for those seeking a multi skilled profile, which is technically essential for hybrid parallel programming. Regarding how to chose the right way, figure Fig. 5 suggests a a basis for a systematic selection.

2) Applications selection: As previously explained, the selection of the applications to be used as case studies is very important in order to best capture the attention of the audience and also create a strong connection with reality. The effect of parallel computing should convincing. Thus, the main criteria for selecting the target applications are: not purely fictive, good illustrative potential, high impact, well parallelizable, not to complicated to implement in parallel. For instance, Bioinformatics topics are particularly exciting in the Brazilian context. Indeed, Brazil has a great biodiversity, however knowledge about the genetic resources and the biotechnological potential of Brazilian species is still very limited. Bioinformatics motivations, challenge and perspectives as they are stated by UnB researchers include:

• Great diversity of organisms in Brazil with biotechnological potential • Major infectious diseases among Brazilian population have been the focus of intensive research • Number of genome and transcriptome projects is increasing very fast due to Brazilian government support • Huge amount of data has been generated and has to be processed and analyzed • Need for active scientific collaborations to provide an efficient computing support to genomic research, with an emphasis on pipeline optimization, information storage and recovering.

We now list some applications extracted from our past and current collaborations:

• One of the main problems Bioinformatics is facing now is how to process genomic data at the rate it is being produced. Once a new biological sequence is discovered, its functional/structural characteristics must be established. For this purpose, newly discovered sequence is compared against the sequences that compose genomic databases, in search of similarities. In order to accelerate the production of results, we could use parallel/distributed variants of the exact methods for pairwise or sequence-profile comparisons, employing many computing units in order to reduce the computation time.

• Flight scheduling is one of the earliest steps in airline operations planning. The schedules quality is crucial for airline profitability and operations reliability. The fleet assignment problem consists in deciding which aircraft types (or fleet) will be assigned to the flights. This problem is routinely considered in airline industry, thus the need of a fast computer program to reach good solutions in a reasonable time. • One of the biggest challenge in Life Science is reading and analyzing the genetic code of the genome. Current high-throughput sequencing machines are able to produce hundreds of millions of short sequences (reads) in a single run and the cost of sequencing has decreased dramatically. However, the task of reconstructing genomes from the large volumes of fragmentary read data still remains an algorithmic challenge that is in practice answered by heuristic and non-satisfactory solutions. Exact methods are computationally costly, thus the need of parallel computing.

• The Amazon has a long history of human settlement. The knowledge that Amazonians have about natural plants and their properties are being lost over time. Currently, there are some attempts to store these information by creating documented catalogs of plants and further analyze them. In this context, number of south America countries that share the Amazonian territory joint together to help Amazonians making their ecosystem safer. The amount of collected information is huge, and process them efficiently is best done by parallel computing. Any connected machine can be accessed remotely through the network. This is the case for most HPC clusters for which such an access is commonly made available for external users. This is also the mean by which HPC resources are shared between different entities. Nevertheless, there are several reasons for planning a local HPC cluster:

1) fulfilling routine computing requests from local entities; 2) allowing internal or private computing activities;

3) serving as a support for parallel computing practices; 4) acting as a resources center; 5) getting or claiming a computing resources independence; 6) fulfilling the need or reaching a milestone of an HPC project; 7) housing a multi-purpose server; to name a few. We now describe the way to target this technical achievement and discuss the main items.

The first technical question, maybe implicit, should be "what is the total peak performance that we want for the supercomputer?", followed by practical concerns like design budget, maintenance budget, lifetime, jobs density and frequency, cooling requirement, to name a few. Other secondary questions might focus on the vendor, specific hardware capabilities, hardware flexibility, robustness, suitability to the expected usage. From the answers of these questions, a picture of the global configuration can be drawn. Now, let move to technical aspect considerations regarding the goal of building a proprietary HPC cluster.

The first and major concern should be the generic compute node, which could be a multicore CPU, a mono-socket manycore CPU, or a multi-socket manycore CPU, coupled or not with GPU card(s). An example of compute node could be a {mono/multi}-socket Intel Broadwell-EP [31]. Figure Fig. 6 displays its main specifications. There pros and cons Fig. 6. Intel Broadwell-EP main characteristics for aggregating to many cores within a single node. The main advantage is the one usually considered for hybrid share/distributed memory computing, which is the reduction of the interconnection network. Having less physical nodes is also an acceptable argument. Another one, very important, is energy, both heat and power consumption. Concerning the cons, the situation is like putting all your eggs within the same basket. In case of (even local) failure, it is likely to disable (and loose) the whole node, which thus yields a severe consequence. Another one is at the OS side, user threads (in a multithreaded program) will collude with standard processes. From the sustained performance viewpoint, a heavy concurrency between too many threads threads within the same node might severely affect the program scalability. This bottleneck mainly results from penalizing memory accesses, which is a crucial concern with NUMA configurations.

Next, since the available compute nodes will have to communicate to each other, hence come questions related to the network like speed (bandwidth and latency), topology, and reliability. Considering a simple Ethernet network will be sufficient and operational.

Storage capacity should be evaluated and extended if necessary. Also for this item, backup mechanisms should be addressed.

For software aspects, the choice will focus on operating system, compilers, programming and runtime libraries, debug tools, monitoring tools, and (batch) schedulers.

For a technically detailed step-by-step explanation and tutorial to build up a standard HPC cluster, we suggest [32].

V. CONCLUSION Designing and implementing a parallel computing curriculum has became a common goal. Indeed, the topic of parallel and distributed computing is now vital because of the emergence and pervasiveness of parallel systems. This topic has several directions and a wide range of applications, thus should be considered a good level of conscience. Managing all related activities within an HPC Center is clearly a good idea, which needs to be considered methodologically. The enthusiasm behind is motivated by the active evolution of HPC systems and the increasing demand of efficient computing solutions for numerous applications. Important facts, advices, schematic views, recommendations and descriptions are reported in this paper in the perspective of creating an active academic HPC Center.

Fig. 1 .

 1 Fig. 1. Computing need standing over computers capability

Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .

 2345 Fig. 2. Overview organization of an academic HPC Center