Xiaodong Liu
email: xiaodong.liu@ec-nantes.fr

Julien Réthoré

Antonius Adrianus Lubrecht

An efficient matrix-free preconditioned conjugate gradient based multigrid method for phase field modeling of fracture in heterogeneous materials from 3D images

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Industrial materials, e.g. composite materials, biological materials, architectured materials, become more and more complex both on components and structures. The complexity of these heterogeneous materials can improve their durability. But it leads to many engineering problems, e.g. failure, which is due to the limited comprehension of these materials. The most common failure mode of engineering materials and structures is fracture, which can sometimes have fearful consequences, e.g. the fracture of a bridge. Researchers pay, therefore, close attention to the prevention of crack induced failure in engineering designs. Numerical modeling of fracture can help researchers to understand and predict failure of cracked structures, e.g. cast iron, composite 10 years later. A parallel version can be found in [START_REF] Ashby | A parallel multigrid preconditioned conjugate gradient algorithm for groundwater flow simulations[END_REF]. Applications to large scale problems were solved by the MGCG in [START_REF] Gilles | Multigrid preconditioned conjugate-gradient method for large-scale wave-front reconstruction[END_REF]. The application of MGCG methods on mixed finite element methods was investigated in [START_REF] Wilson | Multigrid preconditioned conjugate-gradient solver for mixed finite-element method[END_REF]. On the other hand, [START_REF] Bank | Sharp estimates for multigrid rates of convergence with general smoothing and acceleration[END_REF] proposed to use the PCG iterator as a smoother for multigrid algorithms, which is the so-called PCGMG method. [START_REF] Braess | On the combination of the multigrid method and conjugate gradients[END_REF] combined the MGCG method and the PCGMG method to achieve a good efficiency for elasticity problems. The ability to avoid the locking effect of multigrid methods using conjugate gradient solvers was confirmed by [START_REF] Braess | On the combination of the multigrid method and conjugate gradients[END_REF]. We therefore propose to use the PCGMG method to perform the image-based phase field modeling in heterogeneous materials. This choice was made for several reasons: (1), the Jacobi iterator based multigrid method can lead to divergence when solving the mechanical part in phase field models due to the presence of cracks, which leads to huge discontinuities; [START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF], replacing the Jacobi smoother by a PCG iterator does not require a amount of work in terms of implementation based on our previous work; (3), a good parallel performance can therefore be ensured with minimum modifications.

The outline of this paper is the following: Section 2 briefly presents the framework of the staggered phase field modeling. Section 3 introduces the implementation of the proposed preconditioned conjugate based multigrid algorithms by using a finite element discretization. The convergence performance analysis is employed in Section 4. Several improvements are also proposed for problems with large variations in this section. The validation of the overall program and its applications are illustrated in Section 5. The snap-back behavior during the phase field modeling is also presented in this section. An automatic load-control strategy for adjusting the displacement increment is investigated for applications in this section. Conclusions are presented at the end of this paper.

Problem statement

Assuming domain Ω ⊂ R D describes a cracked solid, with its boundary ∂Ω and D the space dimension. Let Γ be the crack in Ω as presented in Figure 1a for a D = 2 case. According to [START_REF] Griffith | Vi. the phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london[END_REF][START_REF] Griffith | The theory of rupture[END_REF] and [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], the energy equilibrium in the cracked domain Ω can be presented as:

E(u, Γ) = E u (u, Γ) + E s (Γ) = Ω W u ( (u)) dΩ + g c X D-1 (1)
where E u (u, Γ) is the elastic energy which can be described as the integral of strain energy density W u ( (u)) in the entire domain. = 1 2 (∇u + ∇u T ) denotes the strain with the displacement u. E s (Γ) denotes the energy required to create the crack in Ω according to the Griffith criterion in [START_REF] Griffith | Vi. the phenomena of rupture and flow in solids, Philosophical transactions of the royal society of london[END_REF][START_REF] Griffith | The theory of rupture[END_REF]. It can be expressed as the product of the Hausdorff surface measure X D-1 and the fracture toughness g c which represents the energy required to create a unit cracked surface for D = 3.

Phase field modeling

However, to solve a fracture problem with a sharp crack (Figure 1a) numerically, is complex. Meanwhile, crack initiation, branching or more generally topological changes of the crack are extremely difficult to model. [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] therefore proposed a smooth function instead of a sharp crack. The formulation to represent the diffuse crack illustrated in Figure 1b, is presented for example in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF]. It reads:

     d -2 c ∆d = 0 in Ω d(x) = 1 on Γ ∇d(x) • n = 0 on ∂Ω (2)
where d(x) describes the states of materials, which is defined as a crack phase field in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF], with x ∈ Ω. ∆ is the Laplacian operator. c represents the actual diffuse crack thickness. Equation ( 2) is the Euler-Lagrange equation associated with the variational problem:

d = Arg {inf Γ d } d ∈ {d|d(x) = 1 on Γ ∀x ∈ Γ} where: Γ d = Ω γ(d, ∇d) dΩ (3)
represents the crack surface functional for a 3D case. γ(d, ∇d) is the crack density function per unit volume, defined in this work by:

γ(d, ∇d) = 1 2 c d 2 + c 2 (∇d) 2 (4)
In the phase field modeling, X D-1 denotes the crack surface functional Γ d . Substituting equation ( 1) and (3), we obtain:

E(u, Γ) = Ω W u ( (u)) dΩ + g c Ω γ(d, ∇d) dΩ = Ω W dΩ (5)
where the free energy W reads:

W = W u ( (u)) + g c 2 c d 2 + g c c 2 (∇d) 2 (6)
To correctly compute W u ( (u)), the unilateral contact formulation presented in [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] is adopted in this work. It reads:

W u ( (u)) = g k (d) 1 2 K [tr( )] 2 + G dev : dev tr( ) ≥ 0 1 2 K [tr( )] 2 + g k (d)G dev : dev tr( ) < 0 ( 7
)
where K is the bulk modulus and G denotes the shear modulus. dev = ε -1 3 tr(ε)I is the deviatoric part of strain tensor. tr( ) is the trace of strain tensor. Ψ + and Ψ -represent the positive and negative part of energy, respectively. It reads:

Ψ + = 1 2 K (tr( )) 2 + G dev : dev tr( ) ≥ 0 G dev : dev tr( ) < 0 (8)
and

Ψ -= 0 tr( ) ≥ 0 1 2 Ktr( ) 2 tr( ) < 0 (9)
The degradation function g k (d) is chosen according to [START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF], which reads:

g k (d) = (1 -k)g(d) + k (10)
where

g(d) = (1 -d) 2
The small parameter k is to maintain the well-posedness of the partially broken system as proposed in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF].

The stress σ can therefore be expressed as:

σ = K d tr(ε)I + 2G d ε - 1 3 tr(ε) (11)
where if

Ψ + ≥ Ψ - K d = [(1 -k)g(d) + k] • K tr( ) ≥ 0 K tr( ) < 0 ( 12
)
and

G d = [(1 -k)g(d) + k] • G (13)
while if Ψ + < Ψ -, d is set to 0 and the initial elastic properties K, G are recovered. The governing equations to determine the displacement field in Ω can be written as:

           d = 0 Ψ + < Ψ - ∇ • σ = 0 in Ω u = U 0 on ∂Ω D σ • n = f ext on ∂Ω N (14)
where U 0 is the prescribed displacement and f ext is the surface forces. The Dirichlet and Neumann boundary conditions are prescribed on ∂Ω D and ∂Ω N , respectively.

According to the second law of thermodynamics, we have:

σ : ˙ -Ẇ ≥ 0 ( 15
)
Substituting σ = ∂W/∂ , we obtain:

- ∂W ∂d ḋ ≥ 0 (16)
As the crack propagation is irreversible, ḋ ≥ 0. We obtain:

- ∂W ∂d = 0 (17)
Considering k ≈ 0, we can therefore obtain the governing equations to compute the phase field d:

     2(1 -d)H -gc c (d -l 2 ∆d) = 0 in Ω d(x) = 1 on Γ ∇d(x) • n = 0 on ∂Ω (18)
where H represents the maximum of elastic energy over history as proposed in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF]. It reads:

H(x, t) = max Ψ + ( (x, t)) -Ψ c (19)
with Ψ c = gc 2 c . It is a threshold that we prescribed according to [START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF]. It means the damage will only appear when H ≥ Ψ c .

The overall staggered phase field modeling can be presented as follows:

Loop on time step t

• Compute displacement u with equation ( 14)

• Calculate H with equation [START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF] • Obtain d with equation [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF] end loop

Implementation

To perform image-based phase field modeling in heterogeneous materials, we use the strategy proposed in our previous papers [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF][START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF] using a finite element based multigrid method. However, we will present a more efficient solver instead of the simple Jacobi iterator proposed in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF][START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF].

Finite element discretization

The procedure of finite element discretization for equation ( 14) was proposed in [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. In this paper, we present a detailed finite element discretization for solving equation [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF]. We still propose to associate one elementary node per image voxel. The chosen element type is 8-node linear hexahedron.

The weak form of equation ( 18) is:

Ω 2H + g c c dd * dΩ - Ω g c c ∆dd * dΩ = Ω 2Hd * dΩ (20)
where d * is the test function. Applying partial integration:

Ω 2H + g c c dd * dΩ + Ω g c c ∇d∇d * dΩ - ∂Ω ∇d • nd * dS = Ω 2Hd * dΩ (21)
with ∂Ω ∇d • nd * dS = 0. Equation ( 21) is therefore:

Ω 2H + g c c dd * dΩ + Ω g c c ∇d∇d * dΩ = Ω 2Hd * dΩ (22)
which can be written in the vector form:

F d L = F d R (23)
with

F d L = Ω 2H + gc c dd * dΩ + Ω g c c ∇d∇d * dΩ F d R = Ω 2Hd * dΩ (24)
Applying finite element discretization:

d ≈ N 1 φ i di
where φ i denotes the shape function. N is the number of unknowns. i denotes the node id. di is thus the approximation of d at node i. Supposing:

d * = φ i
Equation [START_REF] Kettler | A multigrid method and a combined multigrid-conjugate gradient method for elliptic problems with strongly discontinuous coefficients in general domains[END_REF] for elementary node j of one 8-node cubic element becomes:

F d Lj = e 8 g=1 8 i=1 2H + gc c φ i φ j di + g c c ∇φ i ∇φ j di F d Rj = e 8 g=1 2Hφ j ( 25
)
where e is the sum over neighboring elements, g is the sum over Gauss integration points, i is the sum over shape function.

Preconditioned conjugate gradient

In the previous part, we proposed to compute the left and the right hand of equation ( 23) at each elementary node, instead of using typical assembled system matrices, because a matrix-free type solver is applied to reduce memory requirements as presented in [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. However, the Jacobi solver used in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF] can not handle huge discontinuities, e.g. divergence occurs for a spherical inclusion mechanical problem with a material property contrast of 1000 (Figure 5).

In the mechanical part of phase field modeling, huge discontinuities, e.g. contrast larger than 10 5 , can be found due to the presence of cracks. Instead of using the Jacobi solver, we therefore propose to use the preconditioned conjugate gradient (PCG) type relaxation in this work.

We proposed to use a diagonal preconditioner as presented in [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. The idea is to compute the diagonal values of system matrix at each node, e.g. the preconditioner to solve equation [START_REF] Bank | Sharp estimates for multigrid rates of convergence with general smoothing and acceleration[END_REF] at node i reads:

M d i = e 8 g=1 8 i=1 2H + g c c φ i φ i + g c c ∇φ i ∇φ i (26)
With this technique, one does not use any matrix, which is the so-called matrix-free type finite element methods.

The generic PCG solver is presented in the following flow chart:

• Initial guess x with Dirichlet boundary conditions

• Compute preconditioner vector M

• r = F R -F L with r = 0 on Dirichlet boundary conditions

• ite = 1 • while ite ≤ number of relaxations 1. z m i,j,k = r m i,j,k /M m i,j,k
2. update z on ghost points. • end

x represents u for mechanical equations and d for phase field equations. (i, j, k) indicates the node ID. m = 1, 2, 3 for the mechanical part and m = 1 for the phase field equation. r represents the residual which is the difference of the right and the left hand side. p, w and z are temporary vectors. α and β are temporary coefficients required for the PCG method. The ghost points mentioned at the second step of the loop, are created to ensure the presence of neighbors in order to perform the integration over the entire support of the nodes of a given MPI domain. They are indispensable to apply the matrix free type FEM, after the MPI domain decomposition.

Multigrid methods

Besides the PCG solvers, intergrid operators are required in multigrid algorithms. Different from the intergrid operators presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF], specific treatment is required to perform phase field modeling.

Equation [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF] indicates that, H, c and g c on the coarse grid are required to perform multigrid methods. In this work, c is considered as a constant, which has the same value on all grids. For g c , which is assigned to elementary nodes, homogenization techniques are used to obtain coarse grid values as presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF]. Besides c and g c on coarse grids, H is also required on coarse grids. As illustrated in equation ( 19), H is computed at each Gauss integration point. Furthermore, H is a history variable that is to be stored over time steps. However, storing H at each Gauss integration point is too expensive. e.g., for N 3 8-node cubic elements, it is 8 times more expensive to store H at the Gauss integration point than storing a single representative value for each element. To save memory space, we thus propose to store H at each element. It is obtained by a simple average of 8 values at Gauss integration points in each element. The integration of the initial phase field problem is thus approximated but, it will be demonstrated in the examples presented in the last Section that it does not introduce any bias in the overall problem resolution. For the H on coarse grids, we proposed to generate them by a simple average of the corresponding elements from the finest grid. It can be computed recursively as presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF].

With H, c and g c on coarse grids, the PCG relaxations can be executed without any obstacles for equation [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF]. To transform the correction and the residual between grids, standard trilinear prolongation and restriction operators are used.

To solve equation ( 14), equations ( 12) and ( 13) have to be updated at each time step, due to the change of d.

As mentioned above, one uses the strategy of one voxel per elementary node. Initial material properties are therefore assigned to each node. However, in finite element methods, what we use is the material property at the Gauss integration points. One proposes to perform an interpolation with the above mentioned shape function. For example for the shear modulus, it reads:

G g = 8 i=1 φ i G i ( 27
)
where G g denotes the shear modulus at a Gauss integration point and G i is the nodal shear modulus. i is the node id. However, when we apply it to equation ( 13), it becomes a little more complicated due to the presence of d. The d at Gauss integration points must firstly be computed:

d g = 8 i=1 φ i d i Φ + ≥ Φ - d g = 0 Φ + < Φ - (28)
where d g denotes d at Gauss integration points. d i is the nodal value. G is then computed at Gauss integration points by equation [START_REF] Ashby | A parallel multigrid preconditioned conjugate gradient algorithm for groundwater flow simulations[END_REF]. Finally, G d at Gauss integration points, reads:

G g d = [(1 -k)(1 -d g ) 2 + k] • G g (29)
Applying the above strategy on the finest grid is straightforward and mandatory. However, to obtain K d and G d on coarse grids requires several specific treatments.

As demonstrated in equations ( 12) and ( 28), one must first compute tr( ), Φ + and Φ - to define K d and G d at each Gauss integration point. To compute tr( ), Φ + and Φ -, the displacement u is indispensable. However, in phase field modeling, the structure is globally under traction or shearing, i.e. tr( ) ≥ 0. Nevertheless, compression, i.e. tr( ) < 0, can be presented in several small zones, or even several nodes on the finest grid. But, on coarse grids, the displacement has been homogenized, one may not find compression on the coarse girds. Using the displacement on coarse grids to obtain tr( ), Φ + and Φ -, can therefore not guarantee a reasonable representation of the original problem on coarse grids. It can cause slow convergence and even divergence in multigrid algorithms. By consequence, we propose to compute these criteria on Gauss integration points at the finest girds. For the K d and G d on coarse grids, the following strategy is proposed:

• Calculate the value of K, G and d at Gauss integration points on the finest grid.

• Compute tr( ), Φ + and Φ -at these Gauss integration points on the finest grid.

• Obtain K d and G d at these Gauss integration points on the finest grid.

• Calculate the average of K d and G d of all Gauss points in each element. These averages, i.e. K e and G e , represent elementary material property on the finest grid.

• Compute K e and G e in each element on coarse girds with homogenization techniques presented in [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF].

The idea of the above strategy is to use the averaged K d and G d to transfer Gauss integration point material property to element. The coarse gird material property is then computed directly from K e and G e on the finest grid. The best representation of the original problem is achieved. Note that our objective is to solve the problem on the finest grid, the existence of coarse grids is to improve convergence rate. That is why on coarse grid, we do not seek an exact solution of the original problem but making several modifications to achieve the best convergence.

The prolongation and restriction operators used for solving equation ( 14) are tri-linear standard operators.

For details of the multigrid algorithm, please refer to our previous work [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF][START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. Here, we just take some notations. The full multigrid (FMG) cycle and V-Cycle are used in this work. Figure 2 presents a 3-level FMG scheme. One performs ν 0 relaxations on level 1, ν 1 relaxations on each level going up and ν 2 relaxations on each level going down. n cy V-Cycles are used on each level. A trilinear interpolation of the solution of level l is applied to obtain the initial solution of level l + 1. The V-Cycle program for the PCGMG method is illustrated in the 1. Compute the coarse grid material property.

FMG-Cycles

2. Carry out ν 2 relaxations with the PCG solver on the finest level l.

3. Inject the solution and restrict the residual to level l -1.

4. Perform ν 2 relaxations on level l -1.

5. Repeat steps 3,4 until the coarsest grid l = 1.

6. Perform ν 0 relaxations on level 1.

7.

Prolong the correction to the next finer level and perform ν 1 relaxations on this level.

8. Repeat step 7 until the finest level. 9. Loop step 2 -7 until obtaining the required relative error.

10. Output results

Overall program

The above sections introduce each part of the phase field modeling. In this section, the entire program to perform the phase field modeling in heterogeneous materials is presented. Besides of the numerical efficiency, the hybrid MPI/OpenMP is also applied to achieve the best performance both on memory and computational time. The entire program is described in the following flow chart.

• Create MPI topology

• Input images and obtain g c , K and G on the finest grid.

• Calculate g c on all coarse grids • t = 0 and initialization of vectors From the above algorithm, we can find that, in step 2.(c), instead of computing the displacement with the updated d and U e , we propose to always compute the solution u t with the same boundary condition U unit and an updated d. Because equation ( 14) is a static linear elastic problem at each time step, its solution u t is proportional to the prescribed displacement. With this strategy, one can minimize changes in the displacement solution on going to the next time step. It can lead to a faster convergence speed when using an iterative solver.

• while t ≤ number of time step 1. if (t=0)

Performance analysis

In this section, the convergence performance of the proposed method is analyzed. Some improvements are also investigated. The parallel performance is also demonstrated at the end of this section.

Convergence performance of the proposed algorithm

The performance of the PCGMG method is firstly analyzed by solving a spherical inclusion problem in a linear elastic domain Ω as illustrated in Figure 3. The cube size is L. Ω is discretized by 128 3 elements. A 6-level multigrid algorithm is applied. The element size on level l is 2 times larger along each direction than on level l + 1. For this first example, we propose to use simple trilinear restriction and prolongation operators to ensure the transfers between grids. The material properties on coarse grids are computed by a simple average of the material property at corresponding nodes (or element, depending on parameters) from the finest grid. The matrix and the inclusion have the same Poisson ratio ν = 0.29. The Young's modulus of the matrix is E M = 233.43 GPa. Supposing the Young's modulus of inclusion E I is smaller than E M , the contrast is defined by:

c = E M E I
The following boundary conditions are prescribed:

           u z = U o on z = L. u z = 0 on z = 0. u x = 0 at node ( L 2 , 0, 0) u y = 0 at node (0, L 2 , 0)
To analyze the convergence performance, the relative residual Rr is computed by the following

Rr = r T • r F T R • F R × 100% ( 30
)
Figure 4 illustrates the convergence performance for the single level PCG solver and the PCGMG method with a material contrast c = 1000. The relative residual to be reached is set to be: Rr ≤ 10 -6 . A 2-1 V-Cycle is used for the PCGMG solver, which means ν 2 = 2 and ν 1 = 1. Supposing the cost of the transform between grids is negligible, the total cost of a 5 V-Cycle PCGMG is equivalent to 17 relaxations on the finest grid. On the other hand, to achieve the required Rr, the single level solver takes 445 relaxations, which is 26 times more expensive than the PCGMG solver. To analyze the locking effect in multigrid methods, the Poisson ratio ν is set to be 0.49 with c = 1000. Figure 6 illustrates the evolution of the relative residual Rr of PCGMG algorithms. The convergence speed is less than the problem with ν = 0.29, but it converges. With the Jacobi based multigrid algorithms, we obtain a divergence from the first relaxation onward.

According to the performance analysis carried out on the spherical inclusion problem, we confirmed that: (1), PCGMG can handle problems with large material property variations which is the case in phase field models; [START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF], PCGMG is able to overcome the locking effects;

(3), PCGMG is much more efficient than the Jacobi solver based multigrid algorithm; (4), The strategy to use a PCGMG algorithm is much more efficient than a single level PCG solver.

However, as presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF], the spherical inclusion is a simple problem with a regular and symmetric geometry. The real microstructure of heterogeneous materials is much more complicated than a spherical inclusion. A standard multigrid method may not be so efficient for handling complex microstructures. A homogenization technique was therefore proposed in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF] to compute material properties on coarse grids. It can guarantee a reasonable representation of the finest grid problem on coarse grids. Here, we use the real tomographic image of nodular graphite cast iron to analyze the convergence performance of such a proposed homogenization technique for solving a problem with a stress singularity.

As presented in Figure 7, a region of interest with 257 3 voxels is taken from the image of graphite cast iron obtained by [START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF]. The voxel size is 5.06 µm. The strategy to obtain the material property from the gray level in the image is presented in [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. The material property of each material illustrated in Table 1 is applied. The convergence performance is analyzed with an initial crack.

In this work, we propose to work in voxel. The size of the entire domain Ω is 256 3 voxels. 256 3

Rr

Convergence performance

Jacobi c=10 Jacobi c=1e3 PCG c=10 PCG c=1e3 PCG c=1e6 PCG c=1e9 PCG c=1e12 which is illustrated in Figure 8. The following boundary conditions are applied :

          
u z = 1 on z = 256 u z = 0 on z = 0 u x = 0 at node (128, 0, 0) u y = 0 at node (0, 128, 0)

The small parameter k is set to 10 -5 , c = 2h with element size h. A 7-level multigrid algorithm is applied. The element size on level l is 2 times larger than on level l + 1. V-Cycles with ν 2 = 2 and ν 1 = 1 are used to solve equation [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF]. However, ν 2 = 4 and ν 1 = 2 are applied for solving equation ( 14) due to higher complexity. A trilinear restriction and prolongation operator is used for solving both equations. A problem is supposed to be converged when Rr ≤ 10 -6 . The problem is solved by using different strategies for estimating material properties on coarse grids. A first simulation is performed in which material properties on coarse grids are obtained by averaging material properties at corresponding nodes from the finest grid. An other strategy is to use the advanced technique proposed in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF] to recursively compute material properties on coarse grids.

Rr

PCGMG with large Poisson ratio

Figure 9a shows that the advanced technique may not give a significant improvement when solving equation [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro ct imaging[END_REF]. Because the best performance is already achieved with a simple average. But for equation [START_REF] Eid | Multiscale analysis of brittle failure in heterogeneous materials[END_REF], as presented in Figure 9b, the advanced technique can distinctly improve the convergence performance, i.e. the simulation is 2 times cheaper. That is because a simple average of material properties on coarse grid can not guarantee a reasonable representation of the original problem for such a complex case. The result is confirmed with the proposal in [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF].

Parallel performance analysis

A hybrid MPI/OpenMP code is built based on our previous work [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. The detailed performance analysis is already presented in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF]for thermal problems. To be rigorous, a brief graphite cast iron problem with a stress singularity. This problem with a stress singularity can be considered as a time step of phase field modeling: one resolution of the phase field problem plus one resolution of the mechanical problem. The computations are performed on the supercomputer Liger of ICI at Ecole Centrale de Nantes. Each processor of Liger computer consists in 12 cores where OpenMP and MPI can be used. Each node has two sockets or processors. Between nodes, MPI is inevitable. Between sockets, MPI is suggested since OpenMP suffers from poor data access patterns when using two sockets. As a consequence of this computer architecture, the number of OpenMP per MPI task is limited to 12. On Liger computer, each node has 128 GB memory space. With the help of the matrix free strategy, 128 GB is enough for solving on one single core the largest numerical example, i.e. 256 3 elements, presented in this work.

To analyze the parallel performance, the wall time is obtained with different configurations (see in Table 2) for the same simulation. Figure 10 shows the speedup as a function of the number (NB) of cores used on this 256 3 elements problem. Within 100 cores, the speedup is almost optimal. When using 384 cores, an efficiency of about 90% can be found, which means the wall time is reduced to 150 s instead of 51567 s with one core. This efficiency confirms to the result obtained in [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF]. In this work, we do not show the efficiency with 1000 cores. This is because the largest numerical example presented in this work has only 256 3 elements. 384 cores are enough to solve such problem. Using more cores can slow down the speedup. Meanwhile, for an entire phase field modeling with such structure, if we have 1000 time step, it will take less than 42 hours. Note that the first time step is the most expensive one. However, if we do not use parallel

Validation and applications

In this section, the validation of the proposed strategy is presented. Two applications are then presented to show the ability of the proposed strategy.

Validation with three dimensional single notched plate

To validate the proposed strategy, a typical crack propagation in a squared shape part with homogeneous material under a tensile loading is illustrated. The three-dimensional mode-I tension test has already been presented in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF] and [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]. However, only the crack pattern can be found in these two papers, the 3D load-displacement curve is not illustrated, which does not allow us to perform a quantitative comparison. In this work, we therefore take the 2D example parameters from the literature ( [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF][START_REF] Riad | Unified phase field model to simulate both intergranular and transgranular failure in polycrystalline aggregates[END_REF]) and try to find an agreement between the 2D and 3D examples.

The dimension of the part and the initial crack are presented in Figure 11. The Young's modulus is set to be 210 GPa. Its Poisson ratio equals to 0.3. The fracture toughness g c is set to 2.7×10 -3 kN•mm -1 . The domain is discretized by 256 × 64 × 256 elements. The element size is about h = 0.0039 mm. c is supposed to be 0.0078 mm, which is twice the element size. The small parameter k equals to 10 -5 . 7 levels are used for the multigrid algorithm. The element size of level l + 1 is two times smaller than on level l. 192 cores are used to perform this simulation.

The following boundary conditions are applied:

u = 0 on z = 0 mm u = {0, 0, U t} on z = 1 mm
As mentioned above, the displacement is always computed with the same prescribed displacement U unit but a new d. This U unit equals to 1 mm in this example. The following increment is applied for this case:

δU = 1e -5 U t ≤ 5e -3 mm δU = 1e -6 U t > 5e -3 mm
At each time step, the solution is supposed to be converged when the relative residual Rr ≤ 10 -6 .

The same material parameters are used according to the 2D single edge notched tensile test presented in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF][START_REF] Riad | Unified phase field model to simulate both intergranular and transgranular failure in polycrystalline aggregates[END_REF]. From the above time step displacement increment, we can find that in this example, the threshold for H is not applied, i.e. Φ c = 0, for comparing with the results from literature. Figure 12 illustrates the crack pattern at different loading step. It confirms to the 3D results presented in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF] and [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF]. From this Figure, we find that the computed crack thickness is smaller than the thickness prescribed to the initial crack. This is because d = 1 are prescribed on three layers of nodes to introduce the initial crack in the model. But l is set to 2h, the computed crack thickness is therefore smaller.

To perform a quantitative 2D-3D comparison, we propose to compute the reaction force F per unit width along the Z direction. It means to divide the reaction force by the width of the domain, i.e. 0.25 mm. The load-displacement curve should then be similar to 2D regardless of 3D effects. Figure 13 illustrates the reaction force F on a unit width along Z direction as a function of the prescribed displacement U t. A maximum reaction force F max of 0.66 kN/mm can be found in this case. The result is similar to the 2D load-displacement curve presented in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF] and [START_REF] Riad | Unified phase field model to simulate both intergranular and transgranular failure in polycrystalline aggregates[END_REF].

Crack initiation and propagation in the structure with a spherical inclusion

In this section, the proposed strategy is applied for the modeling of crack initiation and propagation in a structure with a spherical inclusion as presented in Figure 3. The size of the cube-shaped structure is 128 × 128 × 128 mm 3 . Properties of each material are given in Table 3. One supposes that the crack resides only in the matrix. The g c in the inclusion is therefore supposed to be 100 times greater than in the matrix. 126 3 8-node cubic elements are used to discretize the cube with the element size h = 1 mm. The Multigrid method uses 6 levels. At time step t = 0, FMG cycles are performed. However, V-cycles with ν 2 = 2 and ν 1 = 1 are applied for the other time steps. It is because from one time step to the next time step, the solution should not change a lot. A good initial solution is already obtained. V-Cycles are therefore the best choice. c = 2h is applied. 192 cores are used to perform this simulation. The following boundary conditions are prescribed:

           u z = U t on z = 128 mm u z = 0 on z = 0. u x = 0
at node (64, 0, 0) mm u y = 0 at node (0, 64, 0) mm At each time step, the solution is supposed to be converged when the relative residual Rr ≤ 10 -6 . At t = 0, U e = 0.138 mm is obtained.

For the displacement increment δU , in this work, we propose an automatic load control strategy to better manage the crack propagation. As is well-known, the basic hypothesis of the staggered phase field is that: the solution increment shall be sufficiently small and the system shall be stable to guarantee a good approximation of the monolithic problem. It means that d should not change a lot from one time step to the next one. We therefore propose to compute are the nodal values of d over the entire domain at time step t and t + 1, respectively. δd max thus estimates the maximum value of the variation of d over the entire domain. The automatic time step control system is constructed based on δd max . The principle is to increase δU when δd max is too small, and to reduce δU when δd max is too large.

For this spherical inclusion problem, the following guidelines are proposed as a first attempt:

δU t+1 =            2δU t
δd max ≤ 0.01 δU t 0.01 < δd max < 0.1 0.5δU t δd max ≥ 0.1 0.01U e δU t ≥ 0.01U e where δU t denotes the value of δU at time step t, with δU 1 = 0.01U e at t=1. To avoid a too large displacement increment, the maximum value of δU is limited at 0.01U e.

A fully cracked structure as presented in Figure 14b is obtained at time step t = 91. The top fracture pattern confirms the result presented in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF]. Only one crack pattern can be found in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF], that is because they used irregular meshes and a refinement only on the top. It is well-known that strain concentration can be found on the two poles of the sphere. Multicrack is expected with this symmetric structure without local defects. The Z direction external force versus the displacement response is illustrated by the red solid line in Figure 15. The above result seems to be correct. However, analyzing Figure 14a, one finds that the crack continues to propagate even when δU becomes very close to 0, i.e. 10 -15 mm. Meanwhile, with only 91 time steps, the crack initiation and propagation are obtained even with a such small displacement increment. It means that the crack propagation might be unstable and the crack propagates too fast.

As a consequence, the stability of the above system at each time step is suspect. The snapback behavior during the crack propagation proposed in [START_REF] Gutiérrez | Energy release control for numerical simulations of failure in quasi-brittle solids[END_REF] and [START_REF] Verhoosel | A dissipation-based arc-length method for robust simulation of brittle and ductile failure[END_REF] is therefore considered. The idea is to allow a decrease fo the prescribed displacement during crack propagation, i.e. applying negative δU . The following automatic control strategy is then proposed. For δU t > 0

δU t+1 =                      2δU t
δd max ≤ 0.01 δU t 0.01 < δd max ≤ 0.05 0.2δU t 0.05 < δd max ≤ 0.15 0.1δU t δd max > 0.15 0.01U e δU t ≥ 0.01U e -10 -9 δU t ≤ 10 -9

For δU t < 0 δU t ≥ -10 -9

δU t+1 =                      0.
with δU 1 = 0.01U e at t=1. Figure 17 illustrates the crack at different time steps with automatic load control. Here, 421 time steps are required to reach complete failure. The failure starts at the two poles of the sphere. It is consistent with zones where there are stain concentration on linear elastic problem as obtained in [START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. A significant difference can be found at the final failure between strategies with and without snap-back, especially for the top crack. Note that the non-symmetric crack patterns for the top and the bottom crack are consistent with the symmetry break introduced by the boundary conditions detailed above. The blue dashed line in Figure 15 shows a typical snap-back behavior. It is obtained that the maximum value of the external force is almost the same. It also confirms that during the macroscopically elastic part, the structure is stable. However, once exceeding this phase, the structure becomes unstable and the snap-back behavior is therefore obtained. Figure 16 presents the evolution of δd max as a function of the time step t. To simplify the comparison, the two curves are plotted in the same X scale. It is obvious that without snap-back the δd max remains large after crack initiation, even the load increment vanishes. However, with snap-back, only one pic in δd max can be found at the end of crack propagation. To be quantitative, the average of δd max over time steps is calculated for these two strategies. 0.20 is obtained for the strategy without snap-back. For the strategy with snapback, one obtains 0.06. It confirms that the automatic load control strategy allowing snap-back can guarantee a minimum variation of d for each time step. It thus legitimates the staggered resolution strategy of the phase field problem.

Crack initiation and propagation in cast iron

In this section, the crack propagation is studied in a real image of cast iron as presented in Figure 7. The geometry and material properties were detailed in Table 1. The finite element discretization strategy and the multigrid parameter setting are the same for the cast iron example presented in Section 4. 384 cores are used to perform this simulation. Here, we work with a voxel size of 5.06 µm. Different from Section 4, here, the phase field modeling is performed without an initial crack. It allows one to analyze the crack initiation and propagation. The following boundary conditions are prescribed: The above mentioned automatic load control strategy is applied. It reads: For δU t > 0 with δU 1 = 0.1U e at t=1. U e = 0.064 is obtained at t = 0. The Z direction external force as a function of the prescribed displacement is presented in Figure 18. A snap-back behavior can be found.

           u z = U t on z = 256 u z = 0 on z = 0 u x =
δU t+1 =                          10δU t δd max < 0.
For δU t < 0 δU t+1 =                          0.
Figure 19 illustrates the crack in cast iron at t=100. Different from the crack propagation in homogeneous materials, the crack path is strongly affected by the graphite nodules. As obtained experimentally under cyclic loading in [START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF], the presence of graphite nodules also strongly affected the crack initiation and propagation. The crack surface is therefore no-longer a plane as presented in Figure 20. Meanwhile, multi-cracking can be found. It confirms the ability to 6 Conclusions and perspectives

Conclusions

In this work, a new and efficient numerical strategy is proposed to solve phase field modeling of fracture problems. The matrix-free algorithm is applied to reduce memory requirement and improve parallel performance. PCGMG algorithms permit one to handle large jumps due to the presence of crack. The automatic load control strategy allows us to observe the snap-back behavior during the crack propagation. Several complicated simulations of fracture problems in heterogeneous materials are demonstrated for the first time. The numerical results confirm several phenomena observed during experiments. The micro-macro interactions are well revealed by these complicated numerical examples. It confirms the effect of the microstructure on the macroscopic behavior. This work opens up a new path for solving large scale 3D fracture problems.

Perspectives

The automatic load control strategy proposed in this work is not perfect. It should be adapted for different problems. A more general strategy could be investigated in the future. However, the proposed strategy ensures that the local dissipation does not exceed a given threshold related to the maximum value of the rate of d. In the obtained load-displacement curves, it is observed that our strategy leads sometimes to unrealistic elastic unloading-reloading due to a too high decrease of the displacement increment. But, as there is no spurious energy dissipation during these artifacts, the envelop of the load-displacement curve as well as the damage evolution over time is well controlled and not biased. reference example is taken from [START_REF] Linse | A convergence study of phase-field models for brittle fracture[END_REF]. The geometry of this 1D bar is presented in Figure 21. The same parameter is applied according to [START_REF] Linse | A convergence study of phase-field models for brittle fracture[END_REF]. Its length L equals to 1 mm. The cross section A is set to 1 mm 2 . This bar is composed of a linear-elastic material with E = 10 MPa and g c = 0.1 N/mm. The origin of axis X is at the center of the bar. It holds therefore a domain with [-L/2, L/2]. At x = -L/2, the bar is fixed. A prescribed displacement U t controlled using the proposed strategy is applied at x = L/2. The internal length scale c is set to 0.05L. The element size h equals to c /32.

The following automatic load control system is applied: for δU t > 0

δU t+1 =               
2δU t δd max < 0.001 δU t 0.001 ≤ δd max ≤ 0.00125 0.001δU t δd max > 0.00125 10 -3 δU t > 10 -3 -10 -4 δU t < 10 -5 with δU t = 10 -3 at the first time step. The distribution of d along the bar at d max = 0.99 is presented in Figure 22. The loaddisplacement curve is also obtained (see Figure 23). These two figures show that the result obtained by our strategy is similar to [START_REF] Linse | A convergence study of phase-field models for brittle fracture[END_REF]. A small difference for the maximum load can be marked and also for the shape of d. However, note that different crack density functions are used in this paper and their work.

Figure 1 :

 1 Figure 1: Crack topology approximation in a D = 2 case

3 .

 3 rz new = r T • z 4. β = rz new /rz old with β = 0 when ite = 1 5. p = z + βp 6. w = F L (p) 7. pw = p T • w 8. α = rz new /pw 9. x = x + αp 10. r = r -αw 11. r = 0 on Dirichlet boundary conditions 12. rz old = rz new 13. ite = ite + 1

Figure 2 :

 2 Figure 2: FMG cycles

  (a) Compute material properties, i.e. K e and G e , depend on d on all grids (b) Compute U e where max(W (x)) = Φ c (c) U t = U e 2. else (a) Calculate prescribed displacement U t = U t + δU (b) Compute material properties depend on d on all grids (c) Compute u t always with a prescribed displacement U unit but a new d (d) Compute H with Ut U unit u t (e) Compute the right hand of equation (18) (f) Compute d with the new H (g) output d 3. end if • end

Figure 3 :

 3 Figure 3: Spherical inclusion in a linear elastic domain

Figure 4 :

 4 Figure 4: Convergence of single level PCG and PCGMG algorithms

Figure 5

 5 Figure 5 presents the convergence performance of the Jacobi solver based multigrid and PCG based multigrid for different material property contrasts. From this Figure, we can conclude:(1), A large contrast, e.g. c = 1000 can lead to non-convergence with the Jacobi based multigrid method. (2), PCGMG is more efficient (faster convergence) than Jacobi based multigrid even for small contrast cases, e.g. c = 10. (3), PCGMG is robust even when the contrast is closed to the numerical precision, e.g. c = 10 12 .To analyze the locking effect in multigrid methods, the Poisson ratio ν is set to be 0.49 with c = 1000. Figure6illustrates the evolution of the relative residual Rr of PCGMG algorithms. The convergence speed is less than the problem with ν = 0.29, but it converges. With the Jacobi based multigrid algorithms, we obtain a divergence from the first relaxation onward.According to the performance analysis carried out on the spherical inclusion problem, we confirmed that: (1), PCGMG can handle problems with large material property variations which is the case in phase field models;[START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF], PCGMG is able to overcome the locking effects;(3), PCGMG is much more efficient than the Jacobi solver based multigrid algorithm; (4), The strategy to use a PCGMG algorithm is much more efficient than a single level PCG solver.However, as presented in[START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and[START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF], the spherical inclusion is a simple problem with a regular and symmetric geometry. The real microstructure of heterogeneous materials is much more complicated than a spherical inclusion. A standard multigrid method may not be so efficient for handling complex microstructures. A homogenization technique was therefore proposed in[START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and[START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF] to compute material properties on coarse grids. It can guarantee a reasonable representation of the finest grid problem on coarse grids. Here, we use the real tomographic image of nodular graphite cast iron to analyze the convergence performance of such a proposed homogenization technique for solving a problem with a stress singularity.As presented in Figure7, a region of interest with 257 3 voxels is taken from the image of graphite cast iron obtained by[START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF]. The voxel size is 5.06 µm. The strategy to obtain the material property from the gray level in the image is presented in[START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. The material property of each material illustrated in Table1is applied. The convergence performance is analyzed with an initial crack.In this work, we propose to work in voxel. The size of the entire domain Ω is 256 3 voxels. 256 3 8-node cubic elements are used to discretize Ω. The initial crack is prescribed as:

  Figure 5 presents the convergence performance of the Jacobi solver based multigrid and PCG based multigrid for different material property contrasts. From this Figure, we can conclude:(1), A large contrast, e.g. c = 1000 can lead to non-convergence with the Jacobi based multigrid method. (2), PCGMG is more efficient (faster convergence) than Jacobi based multigrid even for small contrast cases, e.g. c = 10. (3), PCGMG is robust even when the contrast is closed to the numerical precision, e.g. c = 10 12 .To analyze the locking effect in multigrid methods, the Poisson ratio ν is set to be 0.49 with c = 1000. Figure6illustrates the evolution of the relative residual Rr of PCGMG algorithms. The convergence speed is less than the problem with ν = 0.29, but it converges. With the Jacobi based multigrid algorithms, we obtain a divergence from the first relaxation onward.According to the performance analysis carried out on the spherical inclusion problem, we confirmed that: (1), PCGMG can handle problems with large material property variations which is the case in phase field models;[START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF], PCGMG is able to overcome the locking effects;(3), PCGMG is much more efficient than the Jacobi solver based multigrid algorithm; (4), The strategy to use a PCGMG algorithm is much more efficient than a single level PCG solver.However, as presented in[START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and[START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF], the spherical inclusion is a simple problem with a regular and symmetric geometry. The real microstructure of heterogeneous materials is much more complicated than a spherical inclusion. A standard multigrid method may not be so efficient for handling complex microstructures. A homogenization technique was therefore proposed in[START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF] and[START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF] to compute material properties on coarse grids. It can guarantee a reasonable representation of the finest grid problem on coarse grids. Here, we use the real tomographic image of nodular graphite cast iron to analyze the convergence performance of such a proposed homogenization technique for solving a problem with a stress singularity.As presented in Figure7, a region of interest with 257 3 voxels is taken from the image of graphite cast iron obtained by[START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF]. The voxel size is 5.06 µm. The strategy to obtain the material property from the gray level in the image is presented in[START_REF] Liu | An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using ct images[END_REF]. The material property of each material illustrated in Table1is applied. The convergence performance is analyzed with an initial crack.In this work, we propose to work in voxel. The size of the entire domain Ω is 256 3 voxels. 2563 8-node cubic elements are used to discretize Ω. The initial crack is prescribed as: d = 1 when x < 25.6 & z = 127, 128, 129

Figure 5 :

 5 Figure 5: Convergence performance of Jacobi based multigrid and PCG based multigrid for different material property contrasts

Figure 6 :

 6 Figure 6: Convergence performance of PCGMG solver with ν = 0.49

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: The gray level (GL) in image of cast iron

Figure 10 :

 10 Figure 10: The parallel performance of the hybrid MPI/OpenMP strategy

Figure 11 :

 11 Figure 11: Geometry and initial crack

Figure 12 :

 12 Figure 12: Evolution of crack pattern at different loading step: The isosurface of phase field d with a value of 0.99. The left crack pattern is obtained at U t = 0.005591 mm, the middle one is obtained at U t=0.006 mm, the right one is obtained at U t=0.006273 mm.

Figure 13 :

 13 Figure 13: Reaction force in function of prescribed displacement

  δU in function of t (b) Failure at time step t = 91

Figure 14 :

 14 Figure 14: Displacement increment applied and the fully broken structure

Figure 15 :

 15 Figure 15: Snap-back behavior

Figure 17 :

 17 Figure 17: Crack initiation and propagation in spherical inclusion

Figure 18 :

 18 Figure 18: Load-displacement curve in cast iron until time step = 2390

Figure 19 :

 19 Figure 19: Crack initialization in cast iron. Graphite nodules in green and the crack in red

Figure 20 :Figure 21 :

 2021 Figure 20: Crack surface in cast iron at t = 2390

Figure 22 :

 22 Figure 22: Distribution of the phase field parameter d along the bar at d max = 0.99

Figure 23 :

 23 Figure 23: Load -displacement curves for d max = 0.99

Table 1 :

 1 Material properties in cast iron

	Material	E/GPa ν	g c /kN•mm -1
	Iron	210	0.2 1.73 × 10 -3
	Graphite nodules 21	0.3 1.8×10 -4

Table 2 :

 2 Hybrid setting NB of MPI NB of OpenMP Cores used Nodes used

	1	1	1	1
	8	3	24	1
	16	3	48	2
	64	1	64	4
	32	3	96	4
	32	6	192	8
	64	6	384	16

Table 3 :

 3 Material properties in spherical inclusion problem

	Component E/GPa ν	g c /kN•mm -1
	Matrix	233.43 0.29 2.7 × 10 -3
	Inclusion	2334.3 0.29 2.7×10 -1
	the maximum value of the variation of d from time step t to the time step t + 1 over the entire
	domain. It reads:
		δd max = max(d t+1 i	-d t i )
	where d t i and d t+1 i

  Comparison of δd max between two displacement increment strategies

		0.45
		0.4
		0.35
		0.3
	max	0.25
	d	0.2
		0.15
		0.1
		0.05
		0
	Figure 16:
			0001
		2δU t	0.0001 ≤ δd max < 0.01
		δU t	0.01 ≤ δd max ≤ 0.05
		0.2δU t 0.05 < δd max ≤ 0.15
		0.1δU t δd max > 0.15
		0.1U e δU t ≥ 0.1U e
		-10 -9 δU t ≤ 10 -9

  1δU t δd max < 0.0001 0.5δU t 0.0001 ≤ δd max < 0.01

	δU t	0.01 ≤ δd max ≤ 0.05
	5δU t	0.05 < δd max ≤ 0.15
	10δU t	δd max > 0.15
	-0.1U e δU t ≤ 0.1U e
	10 -9	δU t ≥ -10 -9

To validate the proposed load control strategy, a short comparison is proposed using a crack nucleation problem for a one-dimensional bar with a reduced cross section under tension. The

Acknowledgments

The support of Région Pays de la Loire and Nantes Métropole through a grant Connect Talent IDS is gratefully acknowledged.

A A short comparison between the automatic loadcontrol and the arc length strategy

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.


Ali Compaore
email: acompaore64@gmail.com

Roukiatou Nikièma
email: roukiatoun@yahoo.fr

Rasmané Ouédraogo
email: rouedraogo@imf.org

Foreign Aid and Intergenerational Mobility in Africa

Keywords: F35, 055, I24, C35, J62 Foreign Aid, Intergenerational Mobility, Africa

While there is extensive literature examining the growth and development effects of foreign aid, very little attention has been paid to its potential impact on social mobility. Thus, this paper provides the first empirical evidence on the effect of foreign aid on intergenerational educational mobility in Africa. Drawing on a sample of 28 countries over the period 1970-2010 and using the popular and wellknown probit estimator, we find strong evidence that foreign aid raises the likelihood of experiencing upward educational mobility in the region, while the probability of downward educational mobility tends to be lower in countries that receive a high level of foreign aid. These effects mainly operate through the increased financing for education, the improved education system, and policy, as well as improved education conditions. More interestingly, focusing on the sectoral decomposition of total aid receivedi.e., education sector versus the rest of the economy-, the study highlights that foreign aid to the education sector tends to increase the likelihood of upward educational mobility, contrary to aid allocated to the rest of the economy. Our finding suggests that foreign aid has contributed to improving social mobility in African countries.

I. Introduction

At the end of World War II in the 1950s, official development assistance was deemed as pivotal to fight against poverty in the world with the work of authors such as [START_REF] Domar | Capital expansion, rate of growth, and employment[END_REF], Nurkse [START_REF] Makower | Review of Problems of Capital Formation in Underdeveloped Countries. [Review of Review of Problems of Capital Formation in Underdeveloped Countries[END_REF][START_REF] Rosenstein-Rodan | Problems of industrialisation of eastern and south-eastern Europe[END_REF][START_REF] Rosenstein-Rodan | International aid for underdeveloped countries[END_REF]. Indeed, according to these authors, capital accumulation promotes investment which in turn promotes economic growth. Due to insufficient savings, this accumulation of capital in developing countries can be achieved through debt or foreign aid. These analyses justified several transfers of funds as foreign aid to help countries in need. Thus, several theoretical and empirical analyses have attempted to explore the effects of foreign aid on economic growth, on poverty, on inequalities, and on many other aggregates [START_REF] Arndt | Aid, Growth, and Development: Have We Come Full Circle[END_REF][START_REF] Arndt | Assessing foreign aid's long-run contribution to growth and development[END_REF][START_REF] Sachs | The End of Poverty: Economic Possibilities for Our Time Penguin Books[END_REF][START_REF] Stiglitz | Making Globalization Work New York[END_REF]Friedman, 1958;Bauer, 1972;[START_REF] Easterly | Can foreign aid buy growth[END_REF]Easterly, , 2006[START_REF] Easterly | Introduction: Can't take it anymore?[END_REF][START_REF] Moyo | Dead aid: Why aid is not working and how there is a better way for Africa[END_REF][START_REF] Doucouliagos | Aid Effectiveness on Accumulation: A Meta Study[END_REF]Rajan and Subramanian, 2005).

However, despite the very large and extensive existing literature, no study, to the best of our knowledge, has investigated the effects of foreign aid on intergenerational mobility. Indeed, the issue of intergenerational mobility (IGM) or social mobility has been explored in some studies [START_REF] Becker | A theory of intergenerational mobility[END_REF][START_REF] Chetty | The impacts of neighborhoods on intergenerational mobility II: County-level estimates[END_REF]Narayan et al., 2018;Daude and Robana, 2015;[START_REF] Checchi | Intergenerational persistence of educational attainment in Italy[END_REF]Black and Devereux, 2011;[START_REF] Causa | Intergenerational social mobility in OECD countries[END_REF][START_REF] Azam | Like father, like son? Intergenerational educational mobility in India[END_REF][START_REF] Becker | An equilibrium theory of the distribution of income and intergenerational mobility[END_REF]. Intergenerational Mobility is defined by the change in an individual's social position compared to that of his/her parents. This change in social position can be positive (upward mobility) or negative (downward mobility). The possibility for individuals to move up the income ladder, both throughout their life and in relation to his/her parents, plays an important role in the fight against poverty, in the reduction of inequalities, and even in economic growth (World Bank, 2018). As mentioned in [START_REF] Stuhler | A review of intergenerational mobility and its drivers[END_REF], social mobility and inequalities are closely related. Yet income inequality has increased in many countries over the past decades, and research conducted in the 1990s and 2000s showed that these inequalities are much more persistent from one generation to the next. IGM is thus one of the major concerns in developing countries, especially African countries.

Narayan et al. (2018) emphasized that social mobility is much lower, on average, in Africa. It is one of the regions with the lowest social mobility. Although the average IGM has improved across the developing economies since the 1950s, both absolute IGM and relative IGM have stagnated in Africa. In Africa, 35 percent of people born in the 1980s show higher educational mobility than their parents, compared with approximately 60 percent of the same generation in the average economy of East Asia, Latin America, or the Middle East (Narayan et al., 2018).

Indeed, IGM is apprehended in the literature by several indices among which, educational mobility [START_REF] Stuhler | A review of intergenerational mobility and its drivers[END_REF][START_REF] Ouedraogo | Intergenerational Social Mobility in Africa Since 1920[END_REF]. This is explained, on the one hand, by the scarcity and imprecision of data concerning the level of income of individuals in several countries. On the other hand, the literature argues that education is one of the most important channels for social mobility given the close link between education and labor market participation under human capital theory. Thus, the IGM could be approximated by educational mobility. [START_REF] Alesina | Intergenerational Mobility in Africa[END_REF] find great heterogeneity between countries in upward and downward mobility. According to them, the probability that children born to parents with no education will complete primary schooling exceeds 70% in South Africa and Botswana; while this rate is less than 20% for Sudan, Ethiopia, Mozambique, Burkina Faso, and Malawi. Their analysis also reveals substantial differences within countries according to their administrative area (only 5% in Turkana County, in the northwest of the country, and over 85% in Westlands).

Likewise, Azomahou and Yitbarek (2020) also find country heterogeneity meaning that Nigeria, Guinea, Ghana, and Uganda experienced the highest intergenerational mobility in education, and the Comoros and Madagascar the lowest.

The literature on the economic determinants of IGM identifies factors such as capital market failure, credit constraints, economic transformation, etc. (Narayan et al., 2018). Even if these different factors seem to have more or less a link with foreign aid, studies explicitly considering foreign aid as a determinant do not exist in the literature, to our knowledge. Yet, given the dependence of African countries on foreign aid, external aid could potentially influence IGM in these countries. The volume of foreign aid, as well as foreign aid by sector, appear to be an important tool for economic policies on IGM.

This paper aims to analyze the effect of foreign aid on intergenerational mobility. Specifically, it explores the effect of foreign aid on educational mobility on the one hand and the effect of the different channels of aid on educational mobility on the other hand. This paper uses household survey data from the Integrated Public Use Microdata Series (IPUMS) on 28 developing countries in Africa over the period 1970-2010. We estimate a dichotomous choice specification by using Probit models. We explore the channels through which foreign aid affects upward and downward educational mobility. We argue that foreign aid can affect educational achievements through three channels, including providing financing for education, helping to improve the education system and policy, and improving education cond itions.

The results suggest that the higher the level of foreign aid in percentage of GDP, the higher the likelihood of experiencing upward educational mobility. Inversely, the probability of downward educational mobility tends to be lower in countries that receive high levels of foreign aid. The level of foreign aid to the education sector tends to increase (decrease) the likelihood of upward (downward) educational mobility, contrary to foreign aid allocated to the other sectors. Results also suggest that the spending channel holds and therefore foreign aid affects educational mobility by increasing the government capacity to spend on the educational sector.

Also, the effect of foreign aid on educational mobility operates through the channel of social inclusion and equity. Education conditions appear to be a good transmission channel for upward educational mobility, contrary to downward mobility. The findings of the paper imply that foreign aid could be instrumental to move up the education ladder, which is key to reduce income inequalities and poverty in Africa.

The rest of this paper is organized as follows. Section 2 deals with data and methodology.

Section 3 shows the results, while section 4 concludes the paper.

II. Data and Methodology

Methodology

Our objective is to examine the effect of foreign aid received on intergenerational social mobility in African states. 𝑿 𝑖𝑗,𝑡 is a vector of covariates including households' socio-economic and demographic characteristics. More specifically, these covariates include:

• Age: it represents the age of the individual at the time of the survey. We also added age square to capture any potential generational effects;

• Size of family: this variable capture the number of individuals in the household;

• Gender status: it is a binary variable taking the value of 1 if the individual is a female and 0 otherwise;

• Location: it represents the place of living. As the previous variable, this variable also is a 0-1 dummy taking the value of 1 if the individual lives in rural area and 0 otherwise.

• Access to infrastructure: it includes access to electricity and clean water, which are all binary variables taking the value of 1 if the household has access to electricity or clean water and 0 otherwise.

𝑿 𝑖𝑗,𝑡 also includes the level of development proxied by GDP per capita, the quality of institutions and a binary variable taking the value of 1 if the country is in conflict and 0

otherwise. The quality of institution variable is a composite index constructed by taking the simple average of three governance indicators -regulatory quality, the rule of law and corruption control.

To estimate our dichotomous choice specification given in equation ( 1), we rely on the popular and well-known probit estimator. Indeed, this estimator is designed to fits empirical models with a binary dependent variable, assuming that the probability of a positive outcome is determined by the standard normal cumulative distribution function. Thus, we estimate separate probit models to assess the effect of foreign aid on the probability of upward or downward educational mobility in Africa. Table A2 of the appendix summarizes variables' source and definitions.

Data

A. Data Sources

This study focuses on a panel of 28 African economies (i.e., Benin, Botswana, Burkina Faso, Cameroon, Egypt, Ethiopia, Ghana, Guinea, Kenya, Lesotho, Liberia, Malawi, Mali, Mauritius, Morocco, Mozambique, Nigeria, Rwanda, Senegal, Sierra Leone, South Africa, Sudan, South Sudan, Tanzania, Togo, Uganda, and Zambia, Zimbabwe) examined over the period 1970-2010 based on data availability. More precisely, the sample consists of around 23 million parentschildren matched pairs across 28 African countries.

The dataset is compiled from various sources. First, the variable of major interest in our study is foreign aid in percentage of GDP, which is extracted from OECD databases. We also use sectoral aid data from AidData project of William & Mary's Global Research Institute. This database has the specificity and the advantage to quantify and provide reliable and accurate granular aid data from various donors covering an important period. More interestingly, the database provides data on foreign aid commitments by various economic sectors including agriculture, education, infrastructure, water supply and sanitation, etc. allowing specific sectoral aid-effect analyses.

Second, the dependent variable in this study is social mobility across generations in Africa mainly captured here through educational attainment mobility. 4 We rely on the intergenerational educational mobility measure recently developed by [START_REF] Ouedraogo | Intergenerational Social Mobility in Africa Since 1920[END_REF] based on the Integrated Public Use Microdata Series (IPUMS) international dataset -a project hosted at the University of Minnesota. Intergenerational mobility (IM) includes upward and downward mobility both considered in the analysis. For instance, an upward IM refers to the case where a child born or adopted from non-educated parents who has completed primary school, while a child born or adopted from educated parents who has not completed primary school indicates a downward IM. Our dependent variable is a 0-1 dummies equaling one, if the child experienced upward/downward mobility and 0 otherwise. Variables that capture households' socio-economic and demographic characteristics are also extracted from IPUMS database. Furthermore, GDP per capita is from the IMF's World Economic Outlook, while the quality of institution variable is extracted from the World Bank's Worldwide Governance Indicators. Finally, our conflict variable is retrieved form are from the Uppsala Conflict Data

Program (UCDP) provided by the Department of Peace and Conflict Research, Uppsala University.

B. Patterns of Educational Mobility in Africa

Drawing mainly on Ouedraogo and Syrichas (2021) analysis, interesting patterns and important heterogeneities emerge from educational mobility in Africa over the considered period of this study. First, educational mobility in Africa has significantly improved over the period.

Considering 5-year intervals of birth-cohorts, Figure 1 5) measure the IM up, the likelihood that a child aged 14+ born from illiterate parents finishes primary school. Column (2) and (3) shows the IM up for boys and girls respectively, whereas (4) and ( 5) the IM up separated for urban and rural regions. Columns ( 6) to (10) measure the IM down, the likelihood that a child aged 14+ born from literate parents fails to finish primary school. Column ( 7) and ( 8) shows the IM up for boys and girls respectively, whereas ( 9) and (10) the IM down separated for urban and rural regions. The last rows report simple unweighted averages across the 28 countries. Source: See Ouedraogo and Syrichas (2021)

III. Results

The baseline results, the sectoral aid, the transmission channel and the robustness check are presented in this section.

3.1.Baseline Results

We report in Table 2 the baseline results. In columns 1-3, we use the likelihood of upward educational mobility as the dependent variable, while the likelihood of downward mobility is used in columns 4-6. The results show the coefficients associated with foreign aid are positive and strongly significant at the 1 percent level in columns 1-3. This finding suggests that the higher the level of foreign aid in percentage of GDP, the higher the likelihood of experiencing upward educational mobility. Based on column (1), an increase of foreign aid by 1 percent of GDP is associated with an increase in the probability of upward educational mobility by 3.3 percentage points. 5 Inversely, the coefficient associated with foreign aid is negative and significant at the 1 percent level in columns 4-6. This implies that the probability of downward educational mobility tends to be lower in countries that receive a high level of foreign aid. More precisely, an increase of foreign aid by 1 percent of GDP is associated with a decline in the likelihood of downward educational mobility by 1.2 percentage points.

Table 2. Baseline Results: Effects of Foreign Aid on Educational Mobility

Regarding the control variables, the results are broadly in line with expectations. The coefficients associated with family size, female, and rural area are negative and positive in 5 This is calculated as the marginal effect at mean values. The same applies to the other calculated probabilities in the next sections.

(1)

( columns 1-3 and 4-6, respectively. Individuals from big families and rural areas tend to have lower (higher) probability to experience an upward (downward) educational mobility. The reason could be that financial constraints may restrict the ability for parents who have many children to send them to school, while the lack of schooling infrastructure and poverty are impeding educational mobility in rural areas. Being in rural areas is associated with 9.1 % (7.2%) lower (higher) probability to experience upward (downward) mobility. Regarding female gender, the results imply that girls have lower (upper) probability to upgrade (downgrade) in terms of educational achievements compared to men. It is widely documented that women and girls are poorly educated in sub-Saharan African them men [START_REF] Klasen | Does gender inequality reduce growth and development[END_REF][START_REF] Kazandjian | Gender equality and economic diversification[END_REF]. Table 2 also shows that the coefficients associated with access to electricity and water are positive (negative) and significant at the 1 percent level in columns 1-3 (columns 4-6). That said, having access to electricity and water is associated with higher (lower) probability to experience upward (downward) educational mobility, thus underlining the importance of access to basic infrastructure.

In columns ( 2) and ( 5), we control for the quality of institutions, and the coefficient associated with it is positive and negative, respectively. Thus, individuals from countries with good quality of institutions tend to have higher (lower) probability to experience upward (downward) education mobility. As previous studies have shown, good institutions are instrumental for better educational outcomes [START_REF] Meier | Corruption in the education sector: An introduction[END_REF][START_REF] Hallak | Corrupt schools, corrupt universities: What can be done? International Institute for Education Planning Paris[END_REF]. We also include conflict in columns ( 3) and ( 6). The results show that conflict is associated with lower (higher) probability of upward (downward) educational mobility. In countries affected by conflict, it is difficult for children to attend school not only because schools are usually close as confrontations are going on but also the infrastructures are often destroyed, teachers left the neighborhood and public services are absent (IMF, 2019;[START_REF] Yamada | Impacts of long-lasting civil conflicts on education: Evidence from the 2014 Census of Myanmar[END_REF].

Sectoral Aid

We test whether sectoral foreign aid matters. In Table 3, we used total foreign aid in percentage of GDP. However, one may consider that foreign aid allocated to the education sector could have a different effect compared to aid allocated to the other sectors of the economy. To test this assertion, we disaggregated total foreign aid into two categories: aid for the education sector and aid for non-education sectors. The results are reported in Table 3. We find that the disaggregation between education and other sectors matters. In fact, the coefficient associated with foreign aid allocated to the education sector is positive and strongly significant in column

(1) where upward educational mobility is used as the dependent variable, while the coefficient associated with aid to non-education sectors is negative and significant at the 1 percent level.

That said, foreign aid to the education sector tends to increase the likelihood of upward educational mobility, contrary to foreign aid allocated to the other sectors.

Table 3. Results Obtained Using Sectoral Foreign Aid However, the coefficient associated with foreign aid to non-education sectors is very low, suggesting that its impact on upward educational mobility is marginal. An increase of foreign aid to the education sectors by 1 percent of GDP would increase the probability of upward educational mobility by 7.2, while an increase of foreign aid to non-education sectors by the same magnitude of 1 percent of GDP will reduce upward educational mobility by only 1.3 percent. Inversely, the coefficients associated with aid to the education sector and noneducation sectors are negative and positive in column (3), respectively, implying that foreign aid allocated to the education sector tends to reduce the likelihood of downward educational mobility, while foreign aid to non-education sectors is correlated with higher downward educational mobility. Yet, the impact of foreign aid to non-education sectors on downward educational mobility is marginal. An increase of foreign aid to the education sector by 1 percent of GDP can reduce the probability of downward educational mobility by 3 percentage points, while a rise of aid to non-education by 1 percent of GDP will increase the likelihood of downward mobility by only 1.1 percentage points.

Furthermore, we break down foreign aid allocated to the education sector into aid for basic education, aid for secondary and tertiary education[1], and non-specified education level. The latter includes foreign aid aiming to improve the education system as a whole. The results reported in Table 3, columns ( 2) and ( 4), show that foreign aid allocated to any level of education tends to increase (decrease) the likelihood of upward (downward) educational mobility. The coefficients associated with foreign aid for a basic, secondary and tertiary, and non-specified education levels are all positive and significant at the 1 percent level in column

(2) and negative and significant in column (4). However, we find that the coefficient associated with aid for basic education is lower than the one associated with secondary and tertiary education, suggesting foreign aid for the upper education level has a bigger effect than aid for the lower education level. For instance, an increase of foreign aid for basic education by 1 percent of GDP is correlated with a reduction of the likelihood of downward educational mobility by 3.1 percentage points, against a decline of 6.8 percentage points when aid for secondary and tertiary education increases by 1 percent of GDP. In addition, the results show that foreign aid for non-specified education levels has a higher impact on educational mobility than targeted foreign aid to basic, secondary, and tertiary education levels. This finding suggests that providing foreign aid aiming at improving the education system as a whole could yield many more substantial benefits than targeting some specific educational levels.

Transmission Channels

In this section, we explore the channels through which foreign aid affects upward and downward education mobility. We argue that foreign aid can affect educational achievements through three channels, including providing financing for education, helping to improve the education system and policy, and improving education conditions. To investigate whether the three channels hold, we include them in the estimates. Compared to the baseline results, the inclusion of these variables can have three possible effects on the magnitude of the coefficient of foreign aid (leave it unchanged, decrease, or increase). If the coefficient of foreign aid does not change, this suggests that the impact of foreign aid on educational mobility is unrelated to the added variables (and thus they cannot be considered as channels). If the coefficient declines, then the effect of foreign aid on educational mobility operates through the added channels only.

Finally, if the coefficient increases, then foreign aid has a bigger effect on educational mobility than just through the channels; that is, there is also a direct effect of foreign aid on educational mobility.

Regarding the spending channel, one would assume that foreign aid allows governments to increase spending on the education sector, which in turn is instrumental to building or maintaining school infrastructure or training centers, financing for teachers, and associated needs. To test the spending channel, we use total public spending for the education sector. The results are reported in Table 4, which should be compared to Table 2 (column (1) for upward mobility and (4) for downward mobility). We include education spend ing in columns ( 1) and ( 4). We observe that the coefficient associated with foreign aid declines from 0.04 to 0.03 for upward educational mobility, and from -0.005 to -0.002 for downward mobility. This finding suggests that the spending channel holds and therefore foreign aid affects educational mobility by increasing the government capacity to spend on the educational sector.

As for the channel of education system and policy, we hypothesize that donors could help improve the education system by providing expertise or financing technical assistance necessary to define the education policies of beneficiary countries. Helping countries to create reliable and robust education systems and to plan the expansion of school systems and assess teacher training needs is useful to improve education outcomes. To assess this transmission channel, we use the World Bank's CPIA index on social inclusion and equity. Given that there is no specific index for the education sector, we think that social policy is a good proxy as it includes the education sector and broadly human resources. In Table 4, we include the CPIA index for social inclusion and equity in columns (2) and ( 5). We find that the coefficient associated with foreign aid declines from 0.04 (column 1, Table 4) to 0.02 for upward educational mobility, and is not statistically significant for downward educational mobility.

This result implies that the effect of foreign aid on educational mobility operates through the channel of social inclusion and equity.

Finally, on the education conditions, we argue that donors help improve education conditions by providing direct project financing, which usually does not go to the government budget.

These types of financing are non-fungible and donors monitor and have control over the projects, contributing to infrastructure building and potentially improving education conditions and outcomes. We use the average pupil/teach ratio as an indicator of education conditions. We include this pupil/teach ratio in columns (3) and ( 6). Compared to columns (1) and (4) in Table 4, the results in Table 4 show that the coefficient associated with foreign aid declines for upward education mobility, while increasing for downward mobility (column ( 6)), respectively. Thus, education conditions appear to be a good transmission channel for upward educational mobility, contrary to downward mobility. Third, we split total foreign aid between multilateral aid and bilateral aid. The results are displayed in columns 4-5 and 9-10. We find that our findings do not change regardless of whether we use multilateral aid or bilateral aid. The coefficients associated with both multilateral and bilateral aid are positive and strongly significant in columns ( 4) and ( 5), respectively, while negative in columns ( 9) and (10), respectively. That said, higher multilateral and bilateral aid could increase the probability of upward educational mobility while reducing the likelihood of downward mobility. Furthermore, the coefficients associated with multilateral aid are around two times higher than those of bilateral aid.

IV. Conclusion

Intergenerational mobility (upward and downward) has long occupied the debates on the economic and social level as a way to reduce poverty and inequalities in developing countries.

Foreign aid has also been recognized as a means of financing to help developing countries, particularly in social infrastructure such as education. This study attempted to unveil the potential effect of foreign aid on intergenerational mobility in education. Using survey data from the Integrated Public Use Microdata Series (IPUMS) international dataset, covering 28 developing countries in Africa from 1970-2010, the results suggests that foreign aid raises the likelihood of experiencing upward educational mobility in the region, while the probability of downward educational mobility tends to be lower in countries that receive a high level of foreign aid. These effects mainly operate through the increased financing for education, the improved education system, and policy, as well as improved education conditions. More interestingly, focusing on the sectoral decomposition of total aid receivedi.e., education sector versus the rest of the economy-, the study highlights that foreign aid to the education sector tends to increase the likelihood of upward educational mobility, contrary to aid allocated to the rest of the economy. Our finding suggests that foreign aid has contributed to improving social mobility in African countries. In terms of policy implications, this study highlights the need for donors to channel much more aid to the education sector, which would contribute more to reducing poverty and inequality.

  displays the average upward and downward intergenerational mobility in educational attainment across gender (Panel A & B) and households location (Panel C & D). It emerges that the probability of upward mobility increased over the birth-cohorts, while the likelihood of downward mobility recorded a significant decrease, regardless of the gender and the location of individuals. A deep diveanalysis highlights that the probability of upward mobility reached up to 50% for children born after the millennium (Panel A). In addition, the gender gap has considerably crumbled for first cohorts kids, reflecting the decline in gender inequalities (Panel A & B). Moreover, focus on the rural-urban divide, it came out that, on average, children living in rural regions have a 10% lower probability of upward mobility than their urban regions' peers.

Figure 1 .

 1 Figure 1. Intergenerational Mobility evolution in Africa, Education

Table 1

 1 For instance, countries including Burkina Faso, Ethiopia, and South Sudan report an upward IM's likelihood of less than 15% on average. At the same time, Botswana, Nigeria, South Africa, Mauritius, and Zimbabwe recorded an upward IM above 60%, on average.

presents the educational mobility indices for each country. It shows that less than 41% of children from non-educated parents have completed primary education in Africa. Downward IM is considerable since one-quarter of children born with literate parents do not complete primary education. In addition, intergenerational educational mobility considerably contrasts across African countries on average.

Table 1 .

 1 Average Intergenerational Mobility in Education by Country

Note: The table demonstrates country-level estimates of IM up and IM down in education. Columns (1) to (

Table 4 .

 4 Transmission Channels

	3.4. Robustness Checks

Occupational attainment is also considered as a measure of social mobility. SeeOuedraogo and Syrichas (2020) for further discussion. Throughout the paper, we use social mobility referring to educational mobility.

We check the robustness of the results by employing different specifications. The results of these tests are reported in Table 5.

First, we use foreign aid per capita, which has also been used in some studies. The results are in columns (1) and ( 6). We still find that the coefficient associated with foreign aid per capita is positively associated with upward educational mobility and negatively correlated with downward mobility.

Table 5. Results: Robustness Checks

Second, we use different coverage of total foreign aid by breaking it down into grants and loans.

We find that the coefficients associated with both grants and loans are positive and strongly significant at the 1 percent level in columns (2) and (3), and negative in columns ( 7) and (8). This finding suggests both grants and loans could increase the probability of upward educational mobility while reducing the likelihood of downward mobility. We also find that the coefficient associated with loans is higher than the one associated with grants.

(1)

(