Franz Chouly
email: franz.chouly@univ-fcomte.fr

Alexei Lozinski
email: alexei.lozinski@univ-fcomte.fr

Parareal multi-model numerical zoom for parabolic multiscale problems

We present a time-parallel numerical zoom method for parabolic multiscale problems. A fine finite element solver is defined on a patch where multiscale effects are localized. A coarse finite element solver, based on a coarser timespace discretization, and eventually smoother than the fine one, is defined on the whole domain of the partial differential equation. The coupling between fine and coarse solvers is carried out with an adaptation of the Parareal algorithm. We detail our numerical method and illustrate it with a numerical experiment.

Introduction

For multiscale problems where the multiscale effects are localized, a great variety of numerical methods have been proposed that are based on iterations between two grids: a coarse one on the whole domain, and a fine one on a patch which covers the region endowed with multiscale characteristics. To quote a few examples, let us mention the Chimera method [START_REF] Steger | A chimera grid scheme[END_REF], the method of finite element patches [START_REF] Glowinski | Finite element approximation of multi-scale elliptic problems using patches of elements[END_REF] or the numerical zoom [START_REF] Kamga | Numerical zoom for multiscale problems with an application to nuclear waste disposal[END_REF]. Recently a new method has been proposed to couple different models on the global domain and on the patch, where for instance the coarser model can be a smoothed version of the original one [START_REF] Lozinski | Méthodes numériques et modélisation pour certains problèmes multi-échelles[END_REF][START_REF] Laborde | Numerical zoom for multi-scale and multi-model problems[END_REF]. We show in this note that an extension for parabolic time-dependent problems of this multi-model numerical zoom can be derived using the Parareal framework [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF].

The Parareal algorithm allows time-parallel time-integration of evolution equations [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF][START_REF] Bal | A "parareal" time discretization for non-linear PDE's with application to the pricing of an American put[END_REF][START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF]. This is of interest for long-time simulations of systems with a few degrees of freedom, or where the number of processors is so large that standard (spatial) domain decomposition techniques reach their limit because of bandwidth effects [START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF]. The stability and convergence properties of this algorithm have been studied for parabolic problems in, e.g., [START_REF] Gander | Analysis of the parareal time-parallel time-integration method[END_REF]. For the multi-model numerical zoom described in [START_REF] Lozinski | Méthodes numériques et modélisation pour certains problèmes multi-échelles[END_REF][START_REF] Laborde | Numerical zoom for multi-scale and multi-model problems[END_REF], a parareal time-integration presents attractive features since it allows a time-step for solution of the fine problem on the patch that can be finer than the time-step of the coarser problem on the global domain. Moreover, the finer computations, that require some extra cost, can then be parallelized to some extent. Note there exists recent works where the Parareal framework is adapted to multiscale problems, in the ODE context [START_REF] Legoll | A micro-macro parareal algorithm: application to singularly perturbed ordinary differential equations[END_REF] and in the lattice-Boltzmann context [START_REF] Astorino | Multiscale coupling of finite element and lattice Boltzmann methods for time dependent problems[END_REF].

We first describe the space semi-discretization for a simple parabolic multiscale problem, and then detail our Parareal time-integration method. We illustrate finally the potential interest of this new method with a numerical experiment.

The multiscale problem and its space semi-discretization

Let Ω ⊂ R d be an open bounded polygonal domain (d ≥ 1) and note ∂Ω its boundary. We consider the following diffusion problem:

∂u ∂t -∇ • (K∇u) = f in Ω × (0, T ), u| ∂Ω×(0,T ) = 0, u(•, t = 0) = u 0 , (1)
where u 0 is the initial condition, f is the source term, and where the diffusion matrix K = (K ij ), 1 ≤ i, j ≤ d, satisfies α|ξ| 2 ≤ K(x, t)ξ • ξ ≤ β|ξ| 2 for all ξ ∈ R d and (x, t) ∈ Ω × (0, T ) with some constants 0 < α ≤ β < ∞. We are interested here in the situation where the coefficients K ij are highly oscillating in a small subdomain ω ⊂ Ω. We want to develop a finite element method that provides an approximation for the solution u on a relatively coarse mesh on Ω, certainly too coarse to resolve properly the highly oscillating features inside ω, and then corrects in an iterative manner the numerical solution using another much finer local mesh on ω. We introduce thus a coarse mesh T H on Ω with a corresponding finite element space V H and a fine mesh T h on ω, from which is built a corresponding finite element space V h . Let us suppose for the sake of simplicity that these spaces are based on P k -continuous Lagrange finite elements, with k ≥ 1 (but this is not a restrictive assumption). We further denote V 0 H ⊂ V H and V 0 h ⊂ V h the subspaces of functions that vanish on ∂Ω and ∂ω respectively. We introduce also the finite element space M h on ∂ω as the discrete trace space of V h :

M h = {µ h ∈ C 0 (∂ω) : µ h | E ∈ P k (E) on every boundary edge (face) E of the mesh T h }. (2)
A multiscale semi-discretization in space of problem (1) can be written as:

Find u H : [0, T ] → V 0 H , u h : [0, T ] → V h , and λ h : [0, T ] → M h that satisfy at all time t ∈ (0, T ): Ω ∂u H ∂t v H + Ω K∇u H • ∇v H = Ω\ω f v H - ∂ω λ h v H + ω ∂u H ∂t v H + ω K∇u H • ∇v H , ∀v H ∈ V 0 H , (3)
ω ∂u h ∂t v h + ω K∇u h • ∇v h - ∂ω λ h v h = ω f v h , ∀v h ∈ V h , (4)
∂ω

u h µ h = ∂ω u H µ h , ∀µ h ∈ M h , (5)
u H (•, t = 0) = u 0 H , u h (•, t = 0) = u 0 h , (6)
where u 0 H and u 0 h are some discretizations of the initial condition u 0 (for instance the Lagrange interpolant or the L 2 -projection), and where K is some smooth extension of K inside ω, i.e. K = K on Ω \ ω and K smoother than K on ω. We note that u H is typically non unique, however its restriction on Ω \ ω is unique and provides here an approximation to the exact solution u . Inside ω, the exact solution is approximated by u h while u H is purely fictitious there. Finally, λ h approximates the normal of u on ∂ω.

A Parareal time-discretization for multi-model numerical zoom

We are going now to discretize (3)-( 6) in time. We imagine that some relatively big time-step ∆t is sufficient to discretize the coarse part of the solution u H so that, given u 0 H ∈ V H , we will search for u n H ∈ V H (n = 1, 2, . . . , N ) which are approximations to u H (t n ) at t = t n = n∆t (T = N ∆t). On the contrary some finer discretization in time may be needed for the fine part of the solution u h and we leave the derivatives in time in the corresponding equations. Thus, with the implicit Euler scheme we have the following problem for each time-slab

(t n-1 , t n ), n = 1, 2, . . . , N : Find u n H ∈ V 0 H , u h : [t n-1 , t n ] → V h , and λ h : [t n-1 , t n ] → M h such that Ω u n H -u n-1 H ∆t v H + Ω K∇u n H • ∇v H = Ω\ω f (t n )v H - ∂ω λ h (t n )v H + ω u n H -u n-1 H ∆t v H + ω K∇u n H • ∇v H , ∀v H ∈ V 0 H , (7)
ω ∂u h ∂t v h + ω K∇u h • ∇v h - ∂ω λ h v h = ω f v h , ∀v h ∈ V h for all t ∈ (t n-1 , t n ), ( 8
) ∂ω u h µ h = ∂ω u n H µ h , ∀µ h ∈ M h for all t ∈ (t n-1 , t n ), (9)
u h (•, t = t n-1 ) = U n-1 h . (10)
In equation [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF], the quantity U n-1 h ∈ V h is a known initial condition. For a sequential time-integration method, it would be given by U n-1 h = u h (t - n-1 ), n ≥ 1, and U 0 h = u 0 h . Remark 1 Note that the (weak) boundary condition u h (t) = u n H , t ∈ (t n-1 , t n ), is rather arbitrary in the equation [START_REF] Legoll | A micro-macro parareal algorithm: application to singularly perturbed ordinary differential equations[END_REF]. We could use u n-1 H instead of u n H , for example, without loosing the accuracy a priori, or linear interpolation instead of a piecewise constant interpolation in time.

We propose the following iterative algorithm to compute the initial conditions U n h , n = 1, 2, . . . , N -1 together with the coarse components u n H . We assume that, at each stage k of the algorithm (k ≥ 0), we have already some approximations

(u k,n H ) n=0,...,N ⊂ V H , (U k,n h ) n=0,...,N -1 ⊂ V h , (Λ k,n h ) n=1,...,N ⊂ M h that are compatible in the sense U k,n h | ∂ω = P h,∂ω u k,n H , n = 0, 1, 2, . . . , N -1, (11)
where P h,∂ω is the L 2 (∂ω)-projector onto M h . We proceed as follows to find (u

ω ∂u k+1,n h ∂t v h + ω K∇u k+1,n h • ∇v h - ∂ω λ k+1,n h v h = ω f v h , ∀v h ∈ V h for all t ∈ (t n , t n+1 ), ∂ω u k+1,n h µ h = ∂ω u k,n+1 H µ h , ∀µ h ∈ M h for all t ∈ (t n , t n+1 ), u k+1,n h | t=tn = U k,n
h , for all n = 0, 1, 2, . . . , N -1.

Sequential coarse computations and jump propagation

Set, for n = 1, 2, . . . , N -1, S k,n h := u k+1,n-1 h | t=tn -U k,n h ∈ V h , and S k,0 h = 0. (a) Compute the coarse solution u k+1,n H ∈ V H for n = 1, 2, . . . , N : Ω u k+1,n H -u k+1,n-1 H ∆t v H + Ω K∇u k+1,n H • ∇v H = Ω\ω f (t n )v H - ∂ω Λ k,n h v H + ω u k,n H -u k,n-1 H ∆t v H + ω K∇u k,n H • ∇v H , ∀v H ∈ V 0 H . ( 12
) (b) Compute the corrections ∆U k+1,n h ∈ V h and ∆Λ k+1,n h ∈ M h for n = 1, 2, . . . , N : ω ∆U k+1,n h -∆U k+1,n-1 h ∆t v h + ω K∇(∆U k+1,n h ) • ∇v h - ∂ω ∆Λ k+1,n h v h = 1 ∆t ω S k,n-1 h v h , ∀v h ∈ V h , ( 13
) ∂ω ∆U k+1,n h µ h = ∂ω (u k+1,n H -u k,n H )µ h , ∀µ h ∈ M h . (14)

Correction of the initial values

Update for n = 1, 2, . . . , N :

U k+1,n h = u k+1,n-1 h | t=tn + ∆U k+1,n h , (15)
Λ k+1,n h = λ k+1,n-1 h | t=tn + ∆Λ k+1,n h .
Concerning the initialization of the algorithm, at k = 0, we can solve (1) directly on the coarsest grid (time-step ∆t and finite element space V H ), and with K as a diffusion matrix, so as to obtain u 0,n H . Then U 0,n h and Λ 0,n h are obtained by solving sequentially the fine problem on the patch with the coarsest time-step ∆t (but still with the finite element space V h and the multiscale diffusion matrix K) , so that property [START_REF] Lozinski | Méthodes numériques et modélisation pour certains problèmes multi-échelles[END_REF] is ensured initially, while the computational cost of the initialization remains of the same order than the cost of the sequential steps 2 (a) and (b) in our algorithm.

Upon convergence we have hopefully S ∞,n

h = 0 so that u ∞,n-1 h | t=tn = U ∞,n h
and putting all the u ∞,n h together we obtain the exact solution to ( 7)- [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF]. In practice, we can stop the algorithm after a few correction iterations, or when the norm of the jumps (S k,n h ) is below a given threshold.

Note finally that the proposed algorithm preserves the relationship (11): Proposition 3.1 The discrete solution (U k,n h , u k,n H ) satisfies for all k ≥ 0 and n = 0, 1, 2, . . . , N -1:

U k,n h | ∂ω = P h,∂ω (u k,n H ). (16)
Proof. Assume that (11) holds for a given k ≥ 0, then (16) holds also due to both relationships ( 14) and (15). ✷ Remark 2 The computation of the corrections ∆U k+1,n h and ∆Λ k+1,n h effectuated in equation ( 13) is motivated by the following: if we put together all the contributions u k+1,n h into a single function u k+1 h : [0, T ] → V h , we observe that it satisfies (4) with the extra term -

N -1 n=1 δ(t-t n ) ω S k,n h v h
, where δ(•-t n ) is the Dirac mass at t n . Therefore we proceeded exactly as in the first version of the Parareal method described in [START_REF] Lions | Résolution d'EDP par un schéma en temps "pararéel[END_REF] for an ordinary differential equation. Remark 3 Let us consider the case where the equations for u k+1,n h are discretized by implicit Euler on the same coarse mesh ∆t as those for u H , i.e. u k+1,n-1 h

| t=tn are obtained from U k,n-1 h via ω u k+1,n-1 h | t=tn -U k,n-1 h ∆t v h + ω K(∇u k+1,n-1 h | t=tn ) • ∇v h - ∂ω λk+1,n h v h = ω f v h , ∀v h ∈ V h , ∂ω (u k+1,n-1 h | t=tn )µ h = ∂ω u k,n H µ h , ∀µ h ∈ M h .
Here λk+1,n h is in fact an approximation to λ k+1,n-1 h (t n ). Combining the above equation with ( 12), ( 13), ( 14) and (15) we observe that the algorithm can be rewritten in terms of only u k,n H , U k,n h and Λ k,n h :

Ω u k+1,n H -u k+1,n-1 H ∆t v H + Ω K∇u k+1,n H • ∇v H = Ω\ω f (t n )v H - ∂ω Λ k,n h v H + ω u k,n H -u k,n-1 H ∆t v H + ω K∇u k,n H • ∇v H , ∀v H ∈ V 0 H , ω U k+1,n h -U k+1,n-1 h ∆t v h + ω K∇U k+1,n h • ∇v h - ∂ω Λ k+1,n h v h = ω f (t n )v h , ∀v h ∈ V h , ∂ω u k+1,n h µ h = ∂ω u k+1,n H µ h , ∀µ h ∈ M h .

A numerical experiment

We solve problem (1) with our algorithm, and with d = 2, Ω = [-1 4 ; 3 4 ] 2 , ω = [0; 1 2 ] 2 , T = 0.1. The expression of the diffusion matrix is: K(x, y, t) = (K(x, y)✶ ω (x, y) + 1)I, where K(x, y) = 100H(sin 20πx)H(sin 20πy), H is the Heaviside function, ω = [2/10; 3/10] 2 ⊂ ω and I is the identity matrix in two dimensions. The source term is f (x, y, t) = ✶ ω (x, y)(10 + 20 sin(226πt)).

The initial condition is u 0 = 0. Note that the diffusion matrix, resp. the source term, is highly oscillating in space, resp. in time, only on the patch ω. Finite element discretization is carried out with piecewisecontinuous P 1 Lagrange finite elements, and with structured meshes of sizes H = 4h = 1 8 . For timediscretization, the coarse time-step is ∆t = 0.005 (N = 20), and the fine problem ( 8)-( 9) is discretized also with implicit Euler, and a finer time-step δt = ∆t 10 . Moreover, for the coarse solver, we take K(x, y, t) = I. Numerical experiments are carried out with FreeFEM++ [START_REF] Hecht | Freefem++ documentation[END_REF].

We compare the solution to a reference that is obtained through direct time-discretization of Problem (3)-(6), using the finest time-step δt (5 iterations between coarse and fine solvers are carried out at each time-step). The results are depicted Figure 1, where the solution u c (t n ) = u k,n h (x = 1 4 , y = 1 4 ) at the center of the domain Ω is displayed, as a function of time, for the reference and various Parareal iterations k. We also display the coarse solution (left panel, 'Coarse'), which ignores completely the (spatial and temporal) multiscale effects (recall that K is smooth, and ω ⊂ ω). The initial solution (k=0) is far away from the reference since it is computed with the coarsest time-step ∆t. After the first Parareal correction iteration k = 1 (right panel), the fine temporal multiscale effects are recovered, but there are still important jumps at the interface between the piecewise fine solutions. After the second Parareal correction k = 2, we get very close to the sequential reference and there is no visible evolution of the solution when k is increased. After k = 9 iterations, the algorithm is stopped. Finally, the table below presents the evolution of the jumps (S k,n h ) at each correction iteration k, which appear to converge to 0.

Figure 1 .

 1 Figure 1. Evolution of the solution uc at the center of the domain Ω, for various Parareal iterations k and for reference (sequential) solution.

2

 2 (0,T ;ω) -4.91 -4.36 -5.75 -6.53 -7.16 -7.27 -7.78 -7.94 -8.34 -8.57

Acknowledgements

The first author thanks Florian Lemarié for discussion and comments.