M Biquard

A Rouxel

S Lacroix

H Carfantan

A Monmayrant

H Camon

CSSNET: A LEARNING ALGORITHM FOR THE SEGMENTATION OF COMPRESSED HYPERSPECTRAL IMAGES

Keywords: Compressed hyperspectral imaging, Dual-Disperser CASSI, Compressed images segmentation

The paper presents a semantic segmentation method which is directly applicable to compressed hyperspectral images acquired with a dual-disperser CASSI instrument. It introduses an algorithm based on a shallow neural network that exploits the spectral filtering performed by the optical system and the compressed hyperspectral images measured by the detector. Encouraging results that exploit 50 to 100 less data than the whole hyperspectral datacube on PaviaU and IndianPines datasets are presented.

INTRODUCTION

Classic hyperspectral imagers that produce a whole data cube pose a twofold challenge in environmental observation and monitoring from space, firstly because of the quantity of data to be transmitted to the ground, and secondly because of the resources required for post-processing analyses. Compressed hyperspectral acquisition systems, where each pixel of the acquired images combines several spectral components, can alleviate these issues, firstly by reducing the amount of data to be acquired and transmitted. The interpretation of compressed measurements requires the use of algorithms that exploit various assumptions about the observed scene, such as regularity and redundancy, and an a priori knowledge of the spectral filtering performed by the optical system. The classical approach is to reconstruct the complete hyperspectral cube from the compressed data. State-of-the-art methods such as TwIST [START_REF] Bioucas-Dias | A new TwIST: Two-step iterative shrinkage/thresholding al-gorithms for image restoration[END_REF], GAP-TV [START_REF] Yuan | Generalized alternating projection based total variation minimization for compressive sensing[END_REF], DeSCI [START_REF] Liu | Rank minimization for snapshot compressive imaging[END_REF] and λ-net [START_REF] Miao | lambda-net: Reconstruct hyperspectral images from a snapshot measurement[END_REF] obtain impressive results that rely on theoretical work on "Coded-Aperture Snapshot Spectral Imager" systems (CASSI, [START_REF] Jalali | Snapshot compressed sensing: Performance bounds and algorithms[END_REF][START_REF] Boyer | Compressed sensing with structured sparsity and structured acquisition[END_REF]). Nevertheless, in many cases and particularly in remote sensing, the objective of image acquisition is not to obtain the signal in its raw form, but rather to extract given semantic information from the scene. In such cases, the reconstruction is only an intermediate step in data processing, it can be costly by itself, it still requires post-processing of the reconstructed cube, whereas it can by-passed by directly processing the compressed data.

In this paper, we present a semantic segmentation method which is directly applicable to compressed hyperspectral images acquired with a dual-disperser CASSI instrument. By refraining from reconstructing the raw signal, this work is in line with work done by Davenport [START_REF] Davenport | The smashed filter for compressive classification and target recognition[END_REF][START_REF] Davenport | Signal processing with compressive measurements[END_REF], Arguello [START_REF] David | Algorithm for spectral image target detection from compressive measurements using a joint sparsity model and a CASSI architecture[END_REF][START_REF] Bacca | Coupled deep learning coded aperture design for compressive image classification[END_REF][START_REF] Diaz | Adaptive multisensor acquisition via spatial contextual information for compressive spectral image classification[END_REF] and Arce [START_REF] Zhang | Compressive hyperspectral image classification using a 3d coded convolutional neural network[END_REF]. We propose an algorithm based on a shallow neural network that exploits the spectral filtering performed by the optical system and the compressed hyperspectral images measured by the detector. The next section briefly describes the Dual-Disperser CASSI imager. Section 3 introduces the CSSNet architecture, and section 4 presents some results.

DUAL-DISPERSER CASSI

Our proposal is based on the use of a DD-CASSI hyperspectral imager, originally introduced in [START_REF] Gehm | Single-shot compressive spectral imaging with a dual-disperser architecture[END_REF], on the basis of which we have developped a prototype composed of two symmetrical spectro-imagers on both sides of a micromirror array (DMD) [START_REF] Hemsley | Calibration of programmable spectral imager with dual disperser architecture[END_REF]. The first spectro-imager spatially separates and images the different spectral planes of the scene on the DMD. Each spectral plane is therefore filtered by a different piece of the DMD. The second spectro-imager recombines the filtered spectral planes and re-images them on the camera, cancelling any spectral dispersion (figure 1). The essential property of this system and necessary to our approach is the co-localization: whatever the spatial filtering introduced by the DMD, the spectral components present at a point (x, y) of the camera come only from the spectrum present in the scene at the corresponding point. A consequence of co-localization is a direct access to the panchromatic image on the camera by completely opening the DMD. In the more general case of any filtering pattern on the DMD, we measure at each point of the camera a weighted sum of the different components of the spectrum at this point (equation 1). The weights are directly related to the mask programmed on the DMD and the dispersion of spectral planes on it. It is thus possible from the mask and the dispersion of the system to build a filtering cube. The discretized model of the DD-CASSI imager is expressed as follows:

I r,c = n λ -1 k=0 h r,c,k S r,c,k (1)
with I r,c the intensity measured by the pixel (r, c) of the

CSSNET ARCHITECTURE

Related work

A number of works aim at extracting semantic information from the full hyperspectral cube via neural networks.

Given the lack of labeled cubes, this extraction is done from small patches of the scene. Approaches come to classify them using 1D and 2D convolutions [START_REF] Lee | Going deeper with contextual CNN for hyperspectral image classification[END_REF] [16] [START_REF] Ben Hamida | 3-d deep learning approach for remote sensing image classification[END_REF], or 3D convolutions [18] [19]. Others do semantic segmentation on these patches [START_REF] Pan | DSSNet: A simple dilated semantic segmentation network for hyperspectral imagery classification[END_REF].

Recent works seek, as we do, to segment the scene directly from the compressed images. Approaches inspired by compressed sensing segment the scene based on sparcity assumptions [START_REF] Ramirez | Spectral image classification from optimal coded-aperture compressive measurements[END_REF]. Others rely on Bayesian probabilities to propose iterative segmentation [START_REF] Dunlop-Gray | Experimental demonstration of an adaptive architecture for direct spectral imaging classification[END_REF]. Hao Zhang et al. use neural networks to co-optimize scene segmentation and mask generation [START_REF] Zhang | Joint coded aperture optimization and compressive hyperspectral image classification using 3D coded neural network[END_REF].

Network inputs

The network works with patches of size p. To form the input, we need 3 objects: the filter cube with spatial dimensions of the patch H ∈ [0, 1] p×p×n λ , the associated compressed image I ∈ R p×p and the panchromatic image of the scene P ∈ R p×p . The output is O ∈ R p×p×q with q the number of classes. We estimate S ∈ {0, 1, ..., q-1} p×p with Ŝ = (arg max O i,j ) 0≤i,j≤p-1 .

We want to extract the spectral information of I by removing P . But I is only a percentage of the light flux of P and does not have the same variations. We then define I spec = Ī -P with Ī and P the normalized versions of I and P . The input of the network consists of ( H, Īspec , P ), H and Īspec being the normalized versions of H and I spec .

Network architecture

The CSSNet ("Compressed Semantic Segmentation Network") is divided into two parts. The first one processes the spectral information brought by the compressed image and the filtering cube. The second part mutualizes these data with the spatial information extracted from the panchromatic image.

The network is only made of convolution layers on feature maps of size p × p : a sufficient padding is applied to keep the size of these feature maps.

A detailed diagram of the architecture is given in Figure 2.

Spectral processing. This part includes two blocks of convolutions F 1 and F 2 . The treatment performed is :

H 2 = F 1 ( H), H 3 = F 2 ([H 2 , Īspec ]) (2)
where the operation [., .] represents the concatenation along the channels. H is of dimension p × p × n λ : we consider for the input of F 1 that H is of dimension p × p with n λ channels.

Spatial/spectral processing. The architecture used is inspired by the DSSNet [START_REF] Pan | DSSNet: A simple dilated semantic segmentation network for hyperspectral imagery classification[END_REF] used for the segmentation of the complete hyperspectral cube, to which two convolution layers have been added, and the number of channels modified. This part is decomposed into three convolution blocks F 3 , F 3,sym , and F 4 .

H 4 = F 3 (H 3 ), H 5 = F 3,sym ([H 4 , P ]), O = F 5 (H 5 ) (3)
Optimisation. The loss used is the weighted cross entropy loss to evaluate the difference between Ŝ and S:

L CE = - 1 N N s=1 q-1 k=0 w k 0≤i,j≤p t s i,j,k log(p s i,j,k ) (4)
with N the number of patches, p s i,j,k the softmax probability that pixel (i, j) of patch s belongs to class k, t s i,j,k the class label k for pixel (i, j), and w k the inverse of the median frequency, computed over the whole training dataset, of class k. The number inside each convolution layer represents the number of channels. The inputs of the network are: the patch "panchromatic image" normalized P , the patch "compressed image" normalized Īspec and the patch "associated filtering cube" normalized H. The output of the network is Ŝ which associates to each pixel of the patch a vector of probabilities.

Methodology

Datasets used. The network is tested on two datasets: Indian Pines [START_REF] Baumgardner | 220 Band AVIRIS Hyperspectral Image Data Set[END_REF] and Pavia University [START_REF] Gamba | Pavia University hyperspectral dataset from Grupo de Inteligencia Compututacional[END_REF]. To build the database, these two hyperspectral cubes are run through a DD-CASSI system simulator. The characteristics of the two datasets are presented in Figure 3.

Training and inference separation. Training and inference are performed on compressed images from the same hyperspectral cube. The experiments are done using a well defined separation in which no pixel of the hyperspectral cube involved in the inference is seen during the training phase. An area of 10% is reserved for inference, 80% for training (the rest of the patches are located at the border).

Experimental settings. The patches are of size 7×7. For each dataset, 40 compressed images and one panchromatic image are acquired. 80% of the compressed images are used for training and 20% for inference. For inference, the image is scanned with a sliding window of stride 3. The network is trained on 60 epochs with "early stopping". The learning rate is 0.01, the batch size is 100, the weight decay is 5e-4. The optimizer used is SGD (Stochastic Gradient Descent) with a momentum of 0.9.

Metrics. Two metrics are used: the accuracy (OA) and the average f-score (average of the f-scores per class).

Comparison to existing work. As the implementation details leading to the results presented in [START_REF] Zhang | Joint coded aperture optimization and compressive hyperspectral image classification using 3D coded neural network[END_REF] are not explicit, we only compare ourselves to segmentation networks on full hyperspectral cubes: the 3D CNN of [START_REF] Lee | Going deeper with contextual CNN for hyperspectral image classification[END_REF] implemented by [START_REF] Audebert | Deep learning for classification of hyperspectral data: A comparative review[END_REF] and DSSNet (our own implementation). Each network has been trained 5 times, and the average of the scores is displayed in Figure 4. The underperformance of 3D-CNN on Pavia University, representing the state of the art, is probably due to a problem -whose source we did not find -of hyperparameter setting. For each CSSNet training, the metrics are computed for several DMD masks and the average is taken. Note that with the best performing masks, the accuracy on Pavia University averages 86.5. Compared to the full cube segmentation, there is a loss of performance which is natural because we process 50 and 100 less data volume on Pavia University and Indian Pines. But the results are still acceptable, the details of the errors are shown in Figure 5. On spatial-spectral separation. P is introduced only at the end of the network. A consequent part of this one is dedicated only to the extraction of spectral information. We observe, by introducing P at the same time as Īspec (CSSNet (mix) in Figure 1), a less good generalization of the network to inference. The assumption made is that extracting spectral information from H and Īspec is a more complex task than extracting contours from P . By incorporating P earlier, the network gives it too much importance and "over-learns" the spatial part.

Segmentation results

CONCLUSION

In this paper, we have presented an algorithm for semantic segmentation from compressed hyperspectral data. It is based on a neural network architecture taking into account the specificities of the measurement system. Simulation results indicate lower performances than the state of the art on full hyperspectral cube but the number of acquisitions required to obtain the measurements is drastically reduced by nearly two orders of magnitude. The next steps are to perform iterative segmentation, either with Bayesian fusion or using a RNN type architecture, and perform tests on data sets acquired with our prototype [START_REF] Hemsley | Calibration of programmable spectral imager with dual disperser architecture[END_REF].

Fig. 1 .

 1 Fig. 1. Dual-disperser CASSI architecture.

Fig. 2 .

 2 Fig. 2. Architecture of the CSSNet network. Each blue, yellow, green or red rectangle represents a convolution layer.The number inside each convolution layer represents the number of channels. The inputs of the network are: the patch "panchromatic image" normalized P , the patch "compressed image" normalized Īspec and the patch "associated filtering cube" normalized H. The output of the network is Ŝ which associates to each pixel of the patch a vector of probabilities.

Fig. 3 .

 3 Fig. 3. Used datasets

Fig. 4 .

 4 Fig. 4. Comparison of segmentation algorithm results on PaviaU and Indian Pines. The two shaded lines correspond to the reference algorithms to which we compare ourselves. Unlike our networks, these algorithms take the entire hyperspectral cube as input. The next two rows of the table contain the results of two different versions of our segmentation network. CSSNet (mix) corresponds to a version of CSSNet where the panchromatic image P is introduced together with the compressed image Īspec .

Fig. 5 .

 5 Fig. 5. Segmentation results for DSSNet and CSSNet. Each color represents a class, unlabeled areas are in black. Right: the test areas exploited for the numerical evaluations presented Table 4


Duc Quang Bui

Caroline Japhet
email: japhet@math.univ-paris13.fr

Pascal Omnes
email: pascal.omnes@cea.fr

COUPLING PARAREAL WITH OPTIMIZED SCHWARZ WAVEFORM RELAXATION FOR PARABOLIC PROBLEMS *

Keywords: optimized Schwarz waveform relaxation, Domain decomposition, Parareal in time algorithm, Robin transmission conditions, Convergence rates

We propose and analyze a new parallel paradigm that uses both the time and the space directions. The original approach couples the Parareal algorithm with incomplete optimized Schwarz waveform relaxation (OSWR) iterations. The analysis of this coupled method is presented for a one-dimensional advection-reaction-diffusion equation. We also prove a general convergence result for this method via energy estimates. Numerical results for two-dimensional advection-diffusion problems and for a diffusion equation with strong heterogeneities are presented, to illustrate the performance of the coupled Parareal-OSWR algorithm.

1. Introduction. The Parareal algorithm is a time-parallel method that was proposed by Lions, Maday, and Turinici [START_REF] Lions | A "parareal" in time discretization of PDE's[END_REF] to solve evolutionary problems in parallel. The algorithm is constructed using two solvers: the coarse solver, which is fast but not very accurate, and the fine solver, which is slower but more accurate. The long time interval is divided into smaller windows, the fine solver is performed on each window, using some input initial conditions. The outputs are then corrected by the coarse solver and used as the inputs for the next iteration. In the most simplified convergence analysis, the initial value problem was considered to be an ODE, the fine solver was exact and the coarse solver was the Backward Euler method; in that analysis the outputs were computed exactly and so was the error. More complex analyses of the performance of Parareal were also achieved: for nonexact fine solvers for ODEs [START_REF] Bal | Parareal" time discretization for non-linear PDE's with application to the pricing of an american put[END_REF][START_REF] Baffico | Parallel-in-time moleculardynamics simulations[END_REF], for linear ordinary and partial differential equations [START_REF] Gander | Analysis of the Parareal time-parallel time-integration method[END_REF]; for the nonlinear case [START_REF] Gander | Nonlinear convergence analysis for the Parareal algorithm[END_REF]; for improvement of the parallel efficiency [START_REF] Maday | An adaptive Parareal algorithm[END_REF].

In the construction of Parareal, one may choose suitable coarse and fine solvers to accelerate the process, for example by using an iterative method like, e.g., the Schwarz waveform relaxation (SWR) method [START_REF] Gander | Space-time continuous analysis of waveform relaxation for the heat equation[END_REF]. As Parareal itself is also an iterative method, the overall process would then be composed of outer iterations (Parareal) and inner iterations corresponding to the coarse and fine solvers. In order to save CPU resources, we might think of stopping the inner iterative solvers after a small number of iterations, well before convergence, and hope that the overall convergence could be carried through Parareal iterations. This idea was first raised in [START_REF] Maday | The parareal in time iterative solver: a further direction to parallel implementation[END_REF] for general iterative methods; it was then developed for some iterative methods: for the Jacobi method (see [START_REF] Hernandez | Quelques contributions vers la simulation parallèle de la cinétique neutronique et la prise en compte de données observées en temps réel[END_REF]), for SWR methods (see [START_REF] Guetat | Méthode de parallélisation en temps : application aux méthodes de décomposition de domaine[END_REF][START_REF] Gander | A superlinear convergence estimate for the Parareal Schwarz waveform relaxation algorithm[END_REF]). In this contribution, we are more interested in the latter, which will be recalled next.

SWR algorithms are based on a spatial domain decomposition. The spatial domain is decomposed into overlapping or nonoverlapping subdomains; then the original problem is transformed into multidomain space-time sub-problems. Starting from ini-tial fluxes on the space-time interfaces, each subproblem can then be solved in parallel over the whole time range; then the data are exchanged through the interfaces to create better fluxes. The transmission conditions for the fluxes play an important role in the convergence process (especially for advection-reaction-diffusion problems) and several possibilities can be used, e.g., Dirichlet [START_REF] Giladi | Space-time domain decomposition for parabolic problems[END_REF], Robin [START_REF] Martin | Méthodes de décomposition de domaine de type relaxation d'ondes pour des équations de l'océanographie[END_REF][START_REF]An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF], [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF] or Ventcell [START_REF] Japhet | The Optimized Order 2 method: Application to convectiondiffusion problems[END_REF][START_REF] Halpern | Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems[END_REF]. The last two types of conditions can be optimized [START_REF] Japhet | Optimized Krylov-Ventcell Method. Application to convection-diffusion problems[END_REF][START_REF] Japhet | The best interface conditions for domain decomposition methods: absorbing boundary conditions[END_REF][START_REF] Gander | Optimized Schwarz methods[END_REF], and in that case the corresponding SWR method is called optimized Schwarz waveform relaxation (OSWR). The OSWR method was introduced in [START_REF]Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF], and analyzed and extended to linear advection-reaction-diffusion problems in [START_REF] Martin | Méthodes de décomposition de domaine de type relaxation d'ondes pour des équations de l'océanographie[END_REF][START_REF]An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF]. While in the classical Schwarz method, spatial fluxes are exchanged at each time step [START_REF] Gander | Optimized Schwarz methods[END_REF], the space-time fluxes in the SWR or OSWR methods leave a wide range of choices for the discretization method in the time direction in each sub-problem, which is quite useful in practice (see e.g. [START_REF] Halpern | Optimized Schwarz waveform relaxation and discontinuous Galerkin time stepping for heterogeneous problems[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF][START_REF] Hoang | Space-time domain decomposition for reduced fracture models in mixed formulation[END_REF]).

We now explain the suitability of the SWR method. On the one hand, the analysis of this method, as well as its optimization, can be carried out directly at the (continuous) PDE level, which yields insight which is completely independent from the actual discretization that is used in practice. On the other hand, it can be run in parallel (one processor for each subdomain), and this allows a two-level parallelization process: one level in Parareal iterations, and the other in the SWR (or OSWR) iterations. In addition, as we do not run the solvers until convergence, we shall need to keep additional intermediate outputs from the solvers during Parareal iterations: these outputs are fluxes on the interface. Coupling Parareal with SWR with Dirichlet transmission condition was proposed in [START_REF] Guetat | Méthode de parallélisation en temps : application aux méthodes de décomposition de domaine[END_REF][START_REF] Gander | A superlinear convergence estimate for the Parareal Schwarz waveform relaxation algorithm[END_REF], and in [START_REF] Song | Analysis of two new Parareal algorithms based on the Dirichlet-Neumann/Neumann-Neumann waveform relaxation method for the heat equation[END_REF] for Dirichlet-Neumann/Neumann-Neumann waveform relaxation method.

The purpose of this paper is to propose and analyze a coupled method, called a Parareal-OSWR method [START_REF] Bui | New Space-time Domain Decomposition Algorithms Combined with the Parareal Algorithm[END_REF][START_REF] Bui | Coupling Parareal with Optimized Schwarz Waveform Relaxation for Parabolic Problems[END_REF], that uses Parareal with incomplete OSWR iterations (with optimized transmission conditions) for the fine propagator. We are not aware of any convergence analysis of such method. While OSWR converges much faster than Dirichlet SWR, it also amplifies the difficulties in the convergence analysis, which cannot be performed using the same techniques as for Dirichlet SWR. We provide a convergence analysis, based on energy estimates, for one-dimensional parabolic equations and Robin transmission conditions. Numerical experiments are shown in that case, and also for two-dimensional parabolic equations, with Robin or Ventcel conditions, and in the case where OSWR is accelerated by GMRES (see e.g. [START_REF] Japhet | The best interface conditions for domain decomposition methods: absorbing boundary conditions[END_REF][START_REF] Dolean | An Introduction to Domain Decomposition Methods[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF][START_REF] Ali Hassan | A posteriori stopping criteria for optimized Schwarz domain decomposition algorithms in mixed formulations[END_REF]).

Our paper is organized as follows. In Section 2, we introduce the model problem and state some stability and regularity properties of the solutions of the whole spacetime domain problem and of the sub-space-time-domain problems. In Section 3 we recall the Parareal algorithm. In Section 4 we recall and extend some known results about the OSWR method, which will be used for the analysis of the coupled method. In Section 5, we introduce the coupled Parareal-OSWR algorithm and prove a general convergence result for the method via energy estimates in Section 6. Finally, in Section 8, one-dimensional (1D) and two-dimensional (2D) numerical results, comparing the different methods (OSWR, Parareal, and the coupled Parareal-OSWR algorithm), are discussed. In the 1D case, a comparison is also given with the PSWR method introduced in [START_REF] Song | Analysis of two new Parareal algorithms based on the Dirichlet-Neumann/Neumann-Neumann waveform relaxation method for the heat equation[END_REF].

2. Model problem. Denote Lu := ∂ t u -ν∂ xx u + a∂ x u + bu, Ω = R. For T > 0, we consider the following problem

L(u) = f, in Ω × (0, T ), (2.1a) u(•, 0) = u 0 , in Ω, (2.1b)
where ν, a and b are constants, with ν > 0 and b ≥ 0. The source term f and the initial condition u 0 will be specified in Section 2.2.

2.1. Domain decomposition and notation. We consider a decomposition of the domain Ω into two nonoverlapping subdomains

Ω 1 = (-∞, 0), Ω 2 = (0, +∞),
and introduce the Robin interface operator for i = 1, 2 as follows (see [START_REF]Non-homogeneous Boundary Value Problems and Applications[END_REF][START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF]):

B 1 = ν∂ x - a 2 + p 2 , B 2 = -ν∂ x + a 2 + p 2 . (2.2)
Then, problem (2.1) can be reformulated as the following equivalent multidomain problem [START_REF] Lions | On the Schwarz alternating method. III: A variant for nonoverlapping subdomains[END_REF], with f i = f |Ωi , u i = u |Ωi , and u 0,i = u 0|Ω i , i = 1, 2:

Lu i = f i in Ω i × (0, T ), u i (•, 0) = u 0,i in Ω i , (B i u i )(0, •) = ξ i on (0, T ), i = 1, 2, (2.3)
with ξ i := (B i u j )(0, •) on (0, T ), j = 3 -i, i = 1, 2. (2.4)
In (2.3), the transmision condition involves the parameter p (through B i defined in (2.2)), that is a free parameter chosen such that : a) a Robin subdomain problem of type (2.3) is well-posed, b) it leads to a fast converging algorithm (see Section 4). Both issues will be specified later.

In what follows we will use the notation ξ ξ ξ := (ξ 1 , ξ 2 ) for the Robin data on (0, T ) associated with the solution u of (2.1).

Existence and regularity results

. We introduce the following spaces corresponding to the subdomain problems,

X i = H 1 (Ω i ), i = 1, 2, Y = H 1 4 (0, T ),
and the broken Sobolev space

X := {u ∈ L 2 (Ω) : u |Ωi ∈ X i }, equipped with the norm u X = i u |Ωi 2 Xi 1 2
. With the Robin transmission conditions, we will need more regularity in our analysis, in the anisotropic Sobolev spaces H r,s (Ω × (0, T )) = L 2 (0, T ; H r (Ω)) ∩ H s (0, T ; L 2 (Ω)) defined in [START_REF] Lions | Non-homogeneous Boundary Value Problems and Applications[END_REF].

We recall below some useful regularity properties from [START_REF]Non-homogeneous Boundary Value Problems and Applications[END_REF]. 1 (Ω × (0, T )) and there exists a constant C independent of u 0 and f s.t.

Lemma 2.1. (Regularity of problem (2.1) 1 ) If u 0 ∈ H 1 (Ω) and f ∈ L 2 (0, T ; L 2 (Ω)), problem (2.1) has a unique solution u in H 2,
u H 2,1 (Ω×(0,T )) ≤ C( u 0 H 1 (Ω) + f L 2 (0,T ;L 2 (Ω)) ). (2.5) Lemma 2.2. (Regularity of problem (2.3) 2 ) Let i = 1 or i = 2. If u 0,i ∈ X i , f i ∈ L 2 (0, T ; L 2 (Ω i )) and ξ i ∈ Y, problem (2.
3) has a unique solution u i in H2,1 (Ω i × (0, T )) and there exists a constant C independent of u 0 , f , and ξ i s.t.

u i H 2,1 (Ωi×(0,T )) ≤ C( u 0,i Xi + f i L 2 (0,T ;L 2 (Ωi)) + ξ i Y ). Lemma 2.3. (Trace theorem 3 ) If u ∈ H 2,1 (Ω × (0, T )), then u(•, T ) ∈ X , (B i u)(0, •) ∈ Y, i = 1, 2,
(B i u)(0, .) Y ≤ C u H 2,1 (Ω×(0,T )) , i = 1, 2.
Similar estimates hold by replacing Ω by Ω i , u by u i , and X by X i , for i = 1, 2.

3. Parareal Method. The Parareal method introduced in [START_REF] Lions | A "parareal" in time discretization of PDE's[END_REF] is a numerical method designed to solve evolution problems in parallel. It is based on a decomposition in time of (0, T ) into subintervals : (0, T ) = ∪ N -1 n=0 I n , with

I n = (T n , T n+1 ), for 0 ≤ n ≤ N -1, and 0 = T 0 < T 1 < • • • < T N -1 < T N = T.
Over each such interval generically noted as I := (t 0 , t 1 ), it uses two propagation operators :

• G(I, U 0 ) that provides a rough approximation of u(•, t 1 ), where u is the solution of (2.1), with initial condition u(•, t 0 ) = U 0 . • F(I, U 0 ) that provides a more accurate approximation of u(•, t 1 ). For simplicity, we will consider a regular decomposition of (0, T ), i.e., such that T n+1 -T n = ∆T , ∀n ∈ 0, N -1 . The plain Parareal algorithm [START_REF] Lions | A "parareal" in time discretization of PDE's[END_REF] is as follows.

U k+1 n+1 = G(I n , U k+1 n ) + F(I n , U k n ) -G(I n , U k n ), n = 0, 1, . . . , N -1. (3.1)

end for

We denote by u n the solution of the (sequential) fine propagator at time T n :

u n = F((0, T n ), u 0 ), u n = F(I n-1 , u n-1 ), ∀n ∈ 1, N .
In practice, F will be close to exact, and thus, for the analysis presented here, we suppose that F is the exact propagator, i.e., F((t 0 , t 1 ), ũ0 ) = û(t 1 ) where û is the solution of (2.1a) with initial condition ũ0 at t = t 0 . In particular we can identify u n with u(•, T n ), where u is the solution of (2.1).

4. Optimized Schwarz Waveform Relaxation Method. The OSWR algorithm [START_REF] Gander | Optimal convergence for overlapping and nonoverlapping Schwarz waveform relaxation[END_REF][START_REF]Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF][START_REF]An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF], for solving problem (2.1) is a space-time parallel method based on a domain decomposition in space only Ω 1 = (-∞, 0), Ω 2 = (0, +∞).

Let I = (0, T ), V i = L 2 (Ω i × (0, T )), and I + the set of intervals of R + . The method solves iteratively subproblems on Ω 1 × I and Ω 2 × I, exchanging space-time boundary data5 through the Robin operators B 1 and B 2 (defined in (2.2)), where the parameter p is chosen to optimize the convergence factor of the algorithm.

The method is thus defined using a subproblem solution operator and a transmission operator: for i = 1, 2 :

• the solution operator M i (I, f, u 0,i , ξ i ), i = 1, 2, that maps the available Robin condition ξ i , the initial condition u 0,i and source term f to the subdomain solution u i ,

M i : I + × V i × H 1 (Ω i ) × Y → H 2,1 (Ω i × I) (I, f, u 0,i , ξ i ) → u i , (4.1)
where u i is the solution of the following Robin problem in Ω i × I:

Lu i = f in Ω i × I, u i (•, 0) = u 0,i in Ω i , (B i u i )(0, •) = ξ i on I. i = 1, 2, (4.2)
• the transmission operator B i , i = 1, 2, that maps the available neighbor subdomain solution

u j ∈ H 2,1 (Ω j × I), j = 3 -i, to a new Robin datum ξ i ∈ Y : ξ i = (B i u j )(0, •) on I. Using the definition of M i , problem (2.3)-(2.

4) can be rewritten as

u i = M i (I, f, u 0,i , ξ i ), i = 1, 2, (4.3a)
ξ i = (B i u j )(0, •) on I, j = 3 -i, i = 1, 2. (4.3b)
The OSWR algorithm for solving problem (4.3) (or equivalently (2.1)) is as follows.

Algorithm 4.1 (OSWR)

Choose an initial Robin data ξ ξ ξ 0 = (ξ 0 1 , ξ 0 2 ) on I, for example ξ 0 i = (B i u 0,i )(0, •), on (0, T ), i = 1, 2. for = 1, 2, . . . (OSWR iterations) do 1. Solve the local space-time Robin problems by calculating

u i = M i (I, f, u 0,i , ξ -1 i ), i = 1, 2. (4.4)
2. Update the Robin interface term ξ ξ ξ = (ξ 1 , ξ 2 ), with

ξ i = (B i u j )(0, •) on I, j = 3 -i, i = 1, 2. (4.5)
end for Remark 4.1. By definition of M i and u i in (4.4), the interface condition is

(B i u i )(0, •) = ξ -1 i , i = 1, 2. (4.6)
Then, from (4.5), and using that B i = -B j + p, we obtain, for ≥ 1,

ξ i = -ξ -1 j + 2pu j , j = 3 -i, i = 1, 2. (4.7) Let L ∈ N * .
In what follows, we will denote in compact form

(u L , ξ ξ ξ L ) = OSWR L (I, u 0 , ξ ξ ξ 0 ), (4.8)
where

u L ∈ L 2 (Ω × I) with u L | Ωi = u L i , i = 1, 2, and ξ ξ ξ L = (ξ L 1 , ξ L 2 )
are the output after L iterations of algorithm (4.4)-(4.5) with initial condition u 0 and initial Robin datum ξ ξ ξ 0 on I (the dependency with f is omitted in (4.8), to simplify the notation).

4.1. Stability and convergence. We suppose that f ∈ L 2 (0, T ; L 2 (Ω)). For simplicity, we will use the notation • for the L 2 -norm in Ω or in Ω i , i = 1, 2, and

• I for the (L 2 (I)) 2 -norm. Let (u i , ξ i ), i = 1, 2 be defined by (4.4)-(4.5). For the convergence analysis below, we introduce the following notation for the errors, for i = 1, 2 and ≥ 1: If u 0 ∈ H 1 (Ω), a convergence analysis is done in [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF], we recall it in Theorem 4.2-(i). However, in the context of the coupled Parareal-OSWR method in Section 5, incomplete iterations of the OSWR algorithm are performed at each Parareal iteration. This implies that the new initial condition for the OSWR algorithm, through Parareal iterations, will no longer be in H 1 (Ω), but only in X . Therefore we need the following extended result of Theorem 4.2-(ii) below, that will be used to prove the convergence of the coupled Parareal-OSWR method later. In that case, we suppose that the initial condition of Algorithm 4.1, denoted now by ū0 , verifies ū0 ∈ X , and we introduce the additional notation for the error at time t = 0 : e 0 := ū0 -u(•, 0), and e 0,i := e 0 | Ωi , i = 1, 2.

ζ i := ξ i -ξ i ,
Theorem 4.2. Let L ∈ N * . (i) If u 0 ∈ H 1 (Ω) and ξ ξ ξ 0 ∈ Y 2 , then, Algorithm 4.1 is well defined and we have

L =1 1 2 e (., T ) 2 + i e i 2 L 2 (0,T,H 1 (Ωi)) + 1 2p ζ ζ ζ L 2 I = 1 2p ζ ζ ζ 0 2 I . (4.12)
Hence, Algorithm 4.1 converges for p > 0 in L ∞ (0, T ;

L 2 (Ω 1 )) ∩ L 2 (0, T ; H 1 (Ω 1 )) × L ∞ (0, T ; L 2 (Ω 2 )) ∩ L 2 (0, T ; H 1 (Ω 2 ))
to the solution u of (2.1).

(ii) If the initial condition of Algorithm 4.1 is ū0 ∈ X , and if ξ ξ ξ 0 ∈ Y 2 , then, Algorithm 4.1 is well defined and we have

L =1 1 2 e (., T ) 2 + i e i 2 L 2 (0,T,H 1 (Ωi)) + 1 2p ζ ζ ζ L 2 I = L 2 e 0 2 + 1 2p ζ ζ ζ 0 2 I . (4.13)
Proof. The proof of (i) is done in [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF]Theorem 5.15]. Note that the energy estimate (4.12) is a special case of (4.13) with e 0 = 0.

Proof of (ii). We proceed by induction. We have u 0,i ∈ X i independently of . Moreover, let us suppose that ξ -1 i ∈ Y (this is true for = 1). Then from Lemma 2.2 we have u i ∈ H 2,1 (Ω i × (0, T )). Then from Lemma 2.3, or using (4.7) and the trace theorem, we have ξ i ∈ Y, thus Algorithm 4.1 is well defined.

The proof of the energy estimate (4.13) follows the same steps as in [12, Theorem 5.15] and is detailed in Section 9.

Optimized Robin parameters.

In this section we give the methodology to calculate the Robin parameter p involved in the OSWR method. This parameter is chosen to optimize the convergence factor of the algorithm. Recall that one challenge here is that in the context of the coupled Parareal-OSWR method introduced later, the initial condition of Algorithm 4.1, through Parareal iterations, is in X . Thus we have to extend the calculation of the convergence factor of [START_REF] Gander | Nonlinear convergence analysis for the Parareal algorithm[END_REF], Lemma 5.7 (that was done for ζ 0 i ∈ 0 H 3 4 (0, T )). This can be done by extending

ζ 0 i to a function ζ0 i in H 1 4 (R), from which we obtain ζ i ∈ H 1 4 (R), ≥ 1 (by induction).
Then the convergence factor is obtained by using Fourier transform in time.

By linearity of M i , from (4.3a) and (4.4), the error

e i := u i -u, i = 1, 2, at iteration of the OSWR method, satisfies e i = M i (I, 0, 0, ζ -1 i ), with ζ -1 i = ξ -1 i -ξ i .
Equivalently, e i is solution of the following problem

Le i = 0 in Ω i × I, e i (•, 0) = 0 in Ω i , (B i e i )(0, •) = ζ -1 i on I. i = 1, 2, (4.14)
From (4.3b) and (4.5), we have ζ i = B i e j (0, •), i = 1, 2, and thus the transmission conditions on I in (4.14) also reads

(B i e i )(0, •) = (B i e -1 j )(0, •) on I, j = 3 -i, i = 1, 2. (4.15)
In order to use the Fourier transform in time, we have to extend (4.14)-(4.15) to R.

We proceed by induction on : 6 to obtain a function, denoted by ζ0 i , vanishing on (-∞, 0), and on (T, +∞), for i = 1, 2. Then we can extend equations (4.14) to R, and their solutions, denoted by (ẽ 1 1 , ẽ1

• For = 1, problem (4.14) has an initial Robin datum ζ 0 i ∈ Y that can be extended to H 1 4 (R)
2 ), vanish on (-∞, 0) and coincide with (e 1 1 , e 1 2 ) on (0, T ).

• For ≥ 1, if = 1 then ζ -1 i is defined above; else we define the Robin trace ζ -1 i := (B i ẽ -1 j )(0, •) on {0} × R, that vanishes on (-∞, 0) and coincides with (B i e -1 j )(0, •) on {0} × (0, T ), for i = 1, 2. We assume that ζ -1 i ∈ H 1 4 (R), i = 1, 2.

Let us now prove that ζ

i belongs to H 1 4 (R) for i = 1, 2. The subdomain problems (4.14) are extended on Ω i × R as follows:

Lẽ i = 0 in Ω i × R, ẽ i (•, 0) = 0 in Ω i , (B i ẽ i )(0, •) = ζ -1 i on R, i = 1, 2, (4.16)
and their solution vanish on (-∞, 0) and coincide with (e 1 , e 2 ) on (0, T ). In particular, for ≥ 2, by the definition of ζ -1 i above, (4.15) is extended on {0} × R by

(B i ẽ i )(0, •) = (B i ẽ -1 j )(0, •) on R, j = 3 -i, i = 1, 2. (4.17)
Moreover, by energy estimates, we have ẽ i ∈ L 2 (R; H 1 (Ω i )) and, from the strong form of the equation satisfied by

∂ t ẽ i , we deduce that ∂ t ẽ i ∈ L 2 (R; H -1 (Ω i )), i = 1, 2.
Thus, in what follows, we can use Fourier transform in time. Then, we solve in each subdomain the ordinary differential equation

7 iωê i -ν∂ xx ê i + a∂ x ê i + bê i = 0, i = 1, 2, (4.18)
with the characteristic roots

r -= a - √ d 2ν , r + = a + √ d 2ν , d = a 2 + 4ν(b + iω), (4.19)
where √ d is the complex square-root with positive real part :

√ d = d + a 2 + 4νb 2 + i sign(ω) d -a 2 -4νb 2
,
where d := (a 2 + 4νb) 2 + 16ν 2 ω 2 . Thus, Re(r + ) > 0 and Re(r -) < 0, and the solutions ê i ∈ L 2 (Ω i ), i = 1, 2, are given by8

ê 1 = 2 √ d + p ζ -1 1 (ω)e r + x , ê 2 = 2 √ d + p ζ -1 2 (ω)e r -x , ≥ 1. (4.20)
Then, replacing (4.20) in the transmission conditions (4.17) leads to

∀ ≥ 1, ζ 1 ζ 2 = - √ d + p √ d + p ζ -1 2 ζ -1 1 . (4.21)
Thus, from (4.21) we have, for all ω ∈ R, and for ≥ 1 :

| ζ i (ω)| ≤ | ζ -1 3-i (ω)|, for i = 1, 2. Using the induction hypothesis (i.e. ζ -1 i ∈ H 1 4 (R), for i = 1, 2), the above inequality implies that ζ i ∈ H 1 4 (R), for i = 1, 2. Consequently, we get ζ ζ ζ ∈ (H 1 4 (R)) 2 , ∀ ≥ 1. Setting ζ ζ ζ := ( ζ 1 , ζ 2 ), from (4.21), we have ∀ ≥ 2, ζ ζ ζ = ρ 0 (ω, p) ζ ζ ζ -2 , with ρ 0 (ω, p) := - √ d + p √ d + p 2 . (4.22)
From (4.22), by induction on we obtain,

ζ ζ ζ 2 = (ρ 0 (ω, p)) ζ ζ ζ (0) , ∀ ≥ 1. (4.23)
From (4.23), the convergence factor of the algorithm is ρ 0 (ω, p), defined in (4.22). While we have max ω∈R |ρ 0 (ω, p)| = lim ω→∞ |ρ 0 (ω, p)| = 1, we can use the continuous convergence factor ρ 0 (ω, p) to calculate an efficient Robin parameter for the discrete setting (see, e.g., [START_REF] Japhet | The best interface conditions for domain decomposition methods: absorbing boundary conditions[END_REF][START_REF] Japhet | The Optimized Order 2 method: Application to convectiondiffusion problems[END_REF][START_REF]An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF][START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF]).

Indeed, in numerical computations, the frequency ω is bounded, i.e., we have ω min ≤ ω ≤ ω max , where ω max = π ∆t is the largest discrete frequency supported by the numerical time grid, and ω min = π T is smallest frequency relevant to the global time interval. Defining In practice, the minimization problem (4.24) is solved using the fminsearch function in MATLAB.

5. Coupled Parareal-OSWR method. We set (0, T ) = ∪ N -1 n=0 I n as in Section 3. Then the Parareal-OSWR algorithm is defined using the coarse propagator G of Section 3 and the incomplete9 fine propagator OSWR L defined by (4.8) as follows.

Algorithm 5.1 (Coupled Parareal-OSWR) 1. Choose an initial data (U 0 n ) n∈ 0,N with U 0 0 = u 0 and U 0 n an approximation of u(•, T n ), for example U 0 n := G(I n-1 , U 0 n-1 ), for n = 1, 2, . . . , N . 2. Choose an initial Robin data (ξ ξ ξ 0,0 n ) n∈ 0,N -1 , with ξ ξ ξ 0,0 n := (ξ 0,0 1,n , ξ 0,0 2,n ) on I n , for example ξ 0,0 i,n = (B i U 0 n )(0, •), i = 1, 2. for k=0,1,. . . (Parareal iterations) do 1. On each time interval I n , n = 0, 1, . . . , N -1 : Calculate (u k,L n , ξ ξ ξ k,L n ) = OSWR L (I n , U k n , ξ ξ ξ k,0 n ). (5.1)
2. Set U k+1 0 = u 0 and do Parareal corrections:

U k+1 n+1 = u k,L n (•, T n+1 ) + G(I n , U k+1 n ) -G(I n , U k n ). (5.2)
Update the interface term:

ξ ξ ξ k+1,0 n = ξ ξ ξ k,L n . (5.3) end for Remark 5.1. In (5.1), if L = ∞, then u k,L n ∈ H 1 (Ω)
, for all n ≥ 0, k ≥ 0. However, if L < ∞ and chosen small (e.g. L = 2), then u k,L n ∈ X , for all n ≥ 0, k ≥ 0, but there is very few chance that u k,L n be in H 1 (Ω). Thus, using (5.2) and that u 0 ∈ H 1 (Ω) ⊂ X , we will have U k n ∈ X , for all n ≥ 0, k ≥ 0. 6. Convergence of the Parareal-OSWR algorithm. We will consider the convergence in the L 2 (Ω)-norm. As in Section 4, • will stand for the L 2 -norm in Ω or in Ω i , i = 1, 2, and • I the (L 2 (I)) 2 -norm. Let G be the coarse propagator associated to the source term f = 0. We have the following result Theorem 6.1. We suppose that there exists a constant γ 1 such that

G(I n , U ) ≤ γ 1 U , ∀U ∈ L 2 (Ω), ∀n = 0, 1, . . . , N. (6.1) Then, when k → ∞, U k n converges to u(., T n ) in L 2 (Ω). Moreover, when k → ∞, u k, i,n converges to u| Ωi×[Tn,Tn+1] in L 2 (T n , T n+1 , H 1 (Ω)) for all = 1, 2, . . . , L.
In order to prove Theorem 6.1, we introduce the following notation.

Notation for error estimation. Let u be the solution of problem (2.1), and let (u k, n , ξ ξ ξ k, n ) 1≤ ≤L be the sequence of iterates through the OSWR step (5.1). We define for k = 0, 1, . . ., and n = 0, . . . , N -1, the following:

• E k n := U k n -u n , n = 0, . . . , N , where u n = u(., T n ), • e k, n := u k, n -u, the error in L 2 (Ω × I n ) at each iteration inside step (5.1), with e k, i,n := e k, n | Ωi , i = 1, 2, for = 1, . . . , L, • ζ ζ ζ k, n := ξ ξ ξ k, n -ξ ξ ξ n , for = 1, . . . , L, where ξ ξ ξ n = ((B 1 u)(0, I n ), (B 2 u)(0, I n ))
. Then, by linearity, the algorithm on the error reads as follows. Algorithm 6.1 (Coupled Parareal-OSWR algorithm on the error)

1. Define initial data (E 0 n ) n∈ 0,N with E 0 0 = 0, E 0 n := G(I n-1 , U 0 n-1 ) -u n , where U 0 0 = u 0 , n = 1, 2, . . . , N . 2. Define initial Robin data (ζ ζ ζ 0,0 n ) n∈ 0,N , with ζ ζ ζ 0,0 n := (ζ 0,0 1,n , ζ 0,0 2,n ) on I n , where ζ 0,0 i,n = (B i E 0 n )(0, •), i = 1, 2. for k=0,1,. . . (Parareal iterations) do 1. On each time interval I n , n = 0, 1, . . . , N -1 : Calculate (e k,L n , ζ ζ ζ k,L n ) = OSWR L (I n , E k n , ζ ζ ζ k,0 n ). (6.2)
2. Set E k+1 0 = 0 and do Parareal correction: (see Lemma 6.2), and the terms with E k n are treated using hypothesis (6.1) and proving the new estimate of Lemma 6.3. Lemma 6.2. For all n ∈ 0, N -1 , we have

E k+1 n+1 = e k,L n (•, T n+1 ) + G(I n , E k+1 n -E k n ). ( 6
K k=0 e k,L n (., T n+1 ) 2 ≤ L K k=0 E k n 2 + 1 p ζ ζ ζ 0,0 n 2
In . (6.5)

Proof.

Step (6.2) is the OSWR method on I := I n , with initial condition E k n ∈ X (see Remark 5.1). Thus, from (4.13) with e 0 := E k n , e i := e k, i,n ,

ζ ζ ζ := ζ ζ ζ k, n , we obtain L =1 1 2 e k, n (., T n+1 ) 2 + i e k, i,n 2
L 2 (Tn,Tn+1;H 1 (Ωi)) + 1 2p ζ ζ ζ k,L n 2 In = L 2 E k n 2 + 1 2p ζ ζ ζ k,0 n 2 In .
Moreover, in our algorithm, from (6.4) we have

ζ ζ ζ k+1,0 n = ζ ζ ζ k,L n , thus L =1 1 2 e k, n (., T n+1 ) 2 + i e k, i,n 2
L 2 (Tn,Tn+1;H 1 (Ωi)) + 1 2p ζ ζ ζ k+1,0 n 2 In = L 2 E k n 2 + 1 2p ζ ζ ζ k,0 n 2
In . Summing with respect to k, from 0 to K, we get a telescopic sum on the interface and therefore

K k=0 L =1 1 2 e k, n (., T n+1 ) 2 + i e k, i,n 2
L 2 (Tn,Tn+1;H 1 (Ωi)) + 1 2p ζ ζ ζ K+1,0 n 2 In = L 2 K k=0 E k n 2 + 1 2p ζ ζ ζ 0,0 n 2 In , (6.6)
from which we obtain (6.5).

Lemma 6.3. We suppose that G satisfies (6.1). Then, for all n ∈ 0, N -1 ,

K k=0 E k+1 n+1 2 ≤ 8γ 2 1 K+1 k=0 E k n 2 + 2 K k=0 e k,L n (., T n+1 ) 2 . (6.7)
Proof. Using the triangle inequality in (6.3), and then (6.1), we get

E k+1 n+1 2 ≤ 2 G(I n , E k+1 n -E k n ) 2 + 2 e k,L n (., T n+1 ) 2 , ≤ 2γ 2 1 E k+1 n -E k n 2 + 2 e k,L n (., T n+1 ) 2 , ≤ 4γ 2 1 E k+1 n 2 + E k n 2 + 2 e k,L n (., T n+1 ) 2 .
Then, summing with respect to k, from 0 to K, we have

K k=0 E k+1 n+1 2 ≤ 4γ 2 1 K k=0 E k+1 n 2 + E k n 2 + 2 K k=0 e k,L n (., T n+1 ) 2 ,
from which we deduce (6.7).

With these lemmas, we can now prove Theorem 6.1

Proof of Theorem 6.1. From Lemmas 6.2 and 6.3, we get

K k=0 E k+1 n+1 2 ≤ 8γ 2 1 K+1 k=0 E k n 2 + 2L K k=0 E k n 2 + 2 p ζ ζ ζ 0,0 n 2
In .

Setting

γ 2 = 8γ 2 1 +2L and R n = E 0 n+1 2 + 2 p ζ ζ ζ 0,0 n 2
In , we can rewrite this inequality as

K+1 k=0 E k n+1 2 ≤ γ 2 K+1 k=0 E k n 2 + R n , n ∈ 0, N -1 .
From this inequality, by induction, we obtain

K+1 k=0 E k n+1 2 ≤ γ n+1 2 K+1 k=0 E k 0 2 + n j=0 γ j 2 R n-j , n ∈ 0, N -1 .
Using that E k 0 = 0, ∀k ≥ 0, we finally obtain

K+1 k=0 E k n+1 2 ≤ n j=0 γ j 2 R n-j , n ∈ 0, N -1 . (6.8)
Since the right-hand side of (6.8) does not depend on K, this shows that, for a given n ∈ 0, N , the sum

K k=0 E k n
2 is bounded with respect to K. Hence E k n converges to 0 in the L 2 (Ω)-norm when k → ∞. Moreover, from inequality (6.6), we obtain, for all = 1, 2, . . . , L

K k=0 e k, i,n L 2 (Tn,Tn+1;H 1 (Ωi)) ≤ L 2 K k=0 E k n 2 + 1 2p ζ ζ ζ 0,0 n 2
In .

As the sum in the right-hand side of the above inequality is bounded with respect to K, then the sum in the left is also bounded with respect to K. Hence, e k, i,n tends to 0 in L 2 (T n , T n+1 ; H 1 (Ω i )), i.e, u k, i,n → u i,n in L 2 (T n , T n+1 ; H 1 (Ω i )), for any .

Remark 6.4. Note that Assumption 6.1 is very weak, as we do not require the Lipschitz constant to be strictly lower than 1.10 Remark 6.5. The proof of convergence of the nonoverlapping Parareal-OSWR algorithm in Theorem 6.1 (using energy estimates) is done with U k n ∈ X for all n ≥ 0, k ≥ 0 (see Remark (5.1)). This result is obtained without any correction on U k n so that it is more regular (i.e., in H 1 (Ω)). Thus, in practice we do not need a correction on U k n to have a convergent algorithm. Note that in the case of overlapping subdomains, a correction step will be needed; see [START_REF] Bui | New Space-time Domain Decomposition Algorithms Combined with the Parareal Algorithm[END_REF]. Remark 6.6. In Algorithm 5.1, a possible choice for the initial Robin datum on I n is to take ξ 0,0 i,n = (B i U 0 n )(0), i = 1, 2 (i.e., ξ 0,0 i,n is constant on I n ). A better choice, that improves the convergence of the Parareal-OSWR method, is to define V 0 n as a linear interpolation between U 0 n and U 0 n+1 , and then to take

ξ 0,0 i,n = (B i V 0 n )(0), i = 1, 2.
Thus, in what follows the latter case will be considered.

Numerical results in 1D.

In this section, we consider problem (2.1) with Ω =]0, 1[ and Dirichlet boundary conditions at x = 0 and x = 1. A regular advectiondiffusion example is shown, and then an advection-dominated one (where the Parareal method doesn't perform well). A cell-centered finite volume method is used, with a centered discretization for the advection term (for regular advection-diffusion) or an upwind discretization as in [START_REF] Berthe | Méthodes de décomposition de domaine de type relaxation d'ondes optimisée pour l'équation de convection-diffusion instationnaire discrétisée par volumes finis[END_REF] (in the advection dominated case). For the time discretization, the backward Euler method is used, with time step ∆T for the coarse solver.

For the fine solver with L = ∞, as well as for the coarse solver, we use the OSWR method, with a stopping criterion such that the jump of the optimized transmission conditions, measured in the L 2 -norm on the interface, has been reduced below 10 -13 . Otherwise, we will consider L iterations for the fine solver, with different values of L. A comparison will be done with the PSWR method of [START_REF] Gander | Parareal Schwarz waveform relaxation methods[END_REF], which corresponds to the case L = 1 coupled with an incomplete coarse solver based on only one iteration of the OSWR method.

We take T = 1, a = 1, and b = 0. The domain Ω is decomposed into two subdomains Ω 1 =]0, 1/2[ and Ω 2 =]1/2, 1[ and we choose N = 10 time subintervals. The error is measured in the L ∞ (0, T ; X )-norm. Note that we obtain similar results in the L ∞ (Ω × (0, T ))-norm. The cost of the coarse solver is negligible.

7.1. Regular Advection-Diffusion. We set ν = 1 and choose the right-hand side f and the values of the boundary and initial conditions so that the exact solution is given by u(x, t) = exp(-t) sin(πx). The mesh size and time steps (for the fine solver) are ∆x = 1, 25.10 -3 , ∆t = 1, 69.10 -4 respectively.

In Figure 1, we plot the evolution of the relative error between the Parareal-OSWR (or PSWR) solution and the converged Parareal solution, as a function of the number of Parareal iterations. The case L = ∞ (black circle curve) corresponds to L = 48. The horizontal dashed green line represents the relative scheme error divided by 10, and we consider the number of iterations such that the algorithm error is smaller than this value. This is obtained after 7 iterations for Parareal, and after 13, 19, 39, 77 iterations for Parareal-OSWR with L = 8, L = 4, L = 2, L = 1, respectively. These data are reported in Table 1. We observe that the fastest case is L = 4 with 76 OSWR iterations globally, and it performs similarly as L = 1 or PSWR.

Table 2 shows the gain (in term of fine solver iterations), when the fine solvers are performed in parallel, of Parareal, coupled Parareal-OSWR, or PSWR methods compared to the pure OSWR algorithm (first column). We see that for N = 10 time windows, the gain factor is 1.43 for Parareal and 6. Final gain factor (N/ ) -1.43 6.33 6.33

Table 2 Example in 1D (regular advection-diffusion): Gain factor for Parareal, coupled Parareal-OSWR, and PSWR methods compared to the OSWR solver, with N = 10 time windows 7.2. Advection dominated case. We set ν = 0.001, u 0 (x) = x, ∀x ∈ [0, 1], f = 0, and u(0, t) = 0, u(1, t) = 1, ∀t ∈ (0, T ). The mesh size and time steps (for the fine solver) are ∆t = ∆x = 10 -3 .

In Figure 2, we plot the evolution of the relative error between the Parareal-OSWR (or PSWR) solution and the converged Parareal solution, as a function of the number of Parareal iterations. The classical Parareal algorithm, i.e. the case L = ∞ (black circle curve) does not work at all, since it fully needs ten iterations to converge. In consequence, the final gain factor of the best choice, shown in Table 4, decreases compared to the regular advection-diffusion case, so do all other choices. Figure 2 and Table 3 show that PSWR and Parareal-OSWR with L = 1 are the best choice, and are a good alternative method to classical Parareal, when the latter doesn't work. Indeed, we gain a factor 5.45 for Parareal-OSWR or PSWR. 8. Numerical results in 2D. In this section, we give some numerical illustrations of the performances of the coupled Parareal-OSWR method (Algorithm 5.1) with or without Krylov acceleration (see below), in two space dimensions. In Sections 8.1 and 8.2, these results are shown in the context of the NICEM method [START_REF] Japhet | A new interface cement equilibrated mortar method with Ventcel conditions[END_REF], for the space discretization, that enables the use of more efficient transmission conditions in the OSWR method (i.e., Ventcel conditions) with general domain decomposition and meshes. A regular advection-diffusion case is first considered, with a comparison between Robin and Ventcel conditions. 11 Then, an example with boundary layer and vorticies is given, on conforming and nonconforming meshes, and on long and short time intervals (the latter being the worst case for Parareal). In section 8.3 we show an example in an industrial context [START_REF] Ali Hassan | A posteriori stopping criteria for space-time domain decomposition for the heat equation in mixed formulations[END_REF], using mixed finite elements and Robin transmission conditions. In all cases the backward Euler scheme in time is used, with time step ∆T for the coarse solver.

Solver

The multidomain problem (4.3) can actually be reformulated as an interface problem (see [START_REF] Dolean | An Introduction to Domain Decomposition Methods[END_REF][START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF]) that can be solved by various iterative methods, such as block-Jacobi (which corresponds to the OSWR algorithm) or GMRES. The latter being faster, it will be used and called "OSWRG" in what follows, and the corresponding coupled method is called Parareal-OSWRG. In that case L will designate the number of GMRES iteration, in the fine solver. A comparison between Parareal-OSWR and Parareal-OSWRG will be done in Section 8.1.

For the coarse solver, as well as for the fine solver with L = ∞, we use the OSWRG method, with a stopping criterion such that the jump of the optimized transmission conditions, measured in the L 2 -norm on the interface, has been reduced below 10 -13 . Otherwise, we will consider L iterations for the fine solver, with different values of L. In the results below, the cost of the coarse solver is negligible.

Note that in what follows the error is measured in the L ∞ (0, T ; X )-norm, and that we obtain similar results in the L ∞ (Ω × (0, T ))-norm.

A rotating velocity.

We set Ω =]0, 1[×]0, 1[, T = 21, and consider the 2D problem:

∂ t u + ∇ • (a a au) -ν∆u = f in Ω × (0, T ), (8.1a) u = u D on ∂Ω × (0, T ), (8.1b) u(•, 0) = u 0 in Ω, (8.1c)
with a rotating velocity field a a a = (a x , a y ), where a x = -sin π(y -1 2 ) cos π(x -1 2 ) and a y = cos π(y -1

2 ) sin π(x -1 2 ) ; see Figure 3 In what follows we denote by e the relative scheme error, in L ∞ (0, T ; X )-norm, between the converged Parareal solution and the solution u of problem (8.1), given in (8.2). The term "OSWRG iterations" will designate the iterations of the fine solver. We will give two cases for the diffusion and the space domain decomposition.

Case 1. We take ν = 0.05 and a decomposition of Ω into four subdomains as in Figure 3 in the middle. The number of triangles in the whole domain Ω is 8192 and the mesh size in each subdomain is 0.0221.

In Figure 4 (left), we plot the evolution of the relative error between the Parareal-OSWRG solution (with Ventcel conditions) and the converged Parareal solution, as a function of the number of Parareal iterations, for different values of L. The case L = ∞ (black circle curve) corresponds to L = 10 (average value). The horizontal dashed green line represents 0.1e and we consider the number of iterations such that the algorithm error is smaller than this value. This is obtained after 3 iterations for Parareal or Parareal-OSWRG with L = 8, and after 5, 8, 18 iterations for Parareal-OSWRG with L = 4, L = 2, L = 1, respectively. These data are reported in Table 5 (left). We observe that the fastest case is L = 2 with 16 OSWRG iterations globally.

Table 6 shows the gain (in term of fine solver iterations), when the fine solvers are performed in parallel, of the full Parareal (L = 10) or the coupled Parareal-OSWRG methods compared to the OSWRG algorithm for solving problem (8.1). The first column corresponds to the number of OSWRG iterations to reach 0.1e (where e is the relative scheme error), when the OSWRG method is used to solve (8.1). We observe that for N = 21 time windows, we gain a factor 11.2 for Parareal and a factor 21 for Parareal-OSWRG which is the expected parallel efficiency.

Figure 5, we consider four cases for the fine solver: OSWRG with Ventcel (top left), OSWR with Ventcel (top right), OSWRG with Robin (bottom left), and OSWR with Robin (bottom right). The best choice for the latter is L = 4, and L = 2 else.

The gain factor, when the fine solvers are performed in parallel, compared to the corresponding pure OSWR (or OSWRG) algorithm (with Robin or Ventcel) is as follows: 21 for Parareal-OSWRG with Ventcel (i.e., the expected parallel efficiency), 22.05 for Parareal-OSWRG with Robin, 23.86 for Parareal-OSWR with Ventcel, and 25.66 for Parareal-OSWR with Robin (i.e., the gain factor is a little bit better when using OSWR instead of OSWRG).

In what follows (case 2 and Section 8.2), we will consider the fastest method, that is the Parareal-OSWRG algorithm with Ventcel conditions. Case 2. We take ν = 0.1 and a decomposition of Ω into nine subdomains as in Figure 3 (right). The number of triangles in the whole domain Ω is 18432 and the mesh size in each subdomain is 0.0147.

In Figure 4 (right), the Parareal algorithm (black circle curve) corresponds to L = 24. The horizontal dashed green line represents 0.1e, and is reached after 2 iterations for Parareal, and after 3, 5, 12, 30 iterations for Parareal-OSWRG with L = 8, L = 4, L = 2, L = 1, respectively. These data are shown in Table 5 (right) and we see that the fastest case is L = 4, with a total of 20 OSWRG iterations.

Table 7 shows the gain, when the fine solvers are performed in parallel, of the Parareal or the coupled Parareal-OSWRG methods compared to the pure OSWRG algorithm. In the first column we give the number of OSWRG iterations to reach 0.1e, when it is used as a solver for problem (8.1). We observe that for N = 21 time windows, we gain a factor 10.94 for Parareal and a factor 26.25 for Parareal-OSWRG which represents an efficiency strictly greater than 1, since the number of processors is 21. This is notably better than the expected parallel efficiency.

A boundary layer case with vorticies.

We consider problem (8.1) with f = 0, u 0 (x, y) = 1 -x, u D (x, y) = 1 on {x = 0}, and u D (x, y) = 0 elsewhere, and the following velocity field (see [START_REF] Smolarkiewicz | The multi-dimensional crowley advection scheme[END_REF][START_REF] Martin | Méthodes de décomposition de domaine de type relaxation d'ondes pour des équations de l'océanographie[END_REF]): a a a = 0.32π sin(4πx) sin(4πy), cos(4πx) cos(4πy) ; see Figure 6 In what follows we consider a decomposition of Ω into nine subdomains and consider a uniform mesh (case 1) and then a nonconforming mesh adapted to the physics (case 2). The computed Parareal-OSWRG solution at final time t = T (with the nonconforming mesh) is shown in Figure 6 (right).

Case 1. We consider a uniform mesh as in Figure 3 (right), with a mesh size in each subdomain equal to 0.0147. The time step of the fine solver is ∆t F = 0.0156. In Figure 7 (left), we plot the evolution of the relative error in L ∞ (0, T ; X )-norm, between the Parareal-OSWRG solution and the converged Parareal solution, as a function of the number of Parareal iterations, for different values of L. The Parareal algorithm corresponds to L = 10 and we will compare the other cases of L to this case. We do not know the solution of problem (8.1) and thus the relative scheme error e, however we expect that e is between 10 -1 and 10 -3 , and thus 0.1e between 10 -2 and 10 -4 . Consequently, in Figure 7 (left), the horizontal dashed green lines represent three possible values for 0.1e, and Table 8 shows the number of Parareal and total OSWRG iterations to reach these different values, for L = 1, 2, 4, 8, 10. We observe that the fastest case is L = 2.

We see in Figure 7 two convergence regimes, on the one hand for the case L = ∞ where the iterative method is purely the parareal method, and on the other hand for L = 1 where the method is essentially of the Schwarz type.

More precisely, for finite L (=1, 2, 4, 8) we notice that the first iterations converge at the parareal convergence rate, to take into account the temporal component of the error; then the algorithm takes into account the spatial component of the error with convergence curves that looks like those of the pure Schwarz type iterations (see [START_REF] Gander | Optimized Schwarz methods[END_REF]), since the coupled iterations are essentially Schwarz type iterations, this happening all the faster as L is large.

The OSWRG algorithm used as a solver for problem (8.1) will need 22, 36 and 49 iterations to get a relative error in L ∞ (0, T ; X )-norm (between the iterate solution and the converged OSWRG solution) smaller than 10 -2 , 10 -3 and 10 -4 , respectively. Table 9 shows the gain, when the fine solvers are performed in parallel, of the Parareal or the coupled Parareal-OSWRG methods compared to the OSWRG algorithm. In this table we take the values obtained for e = 10 -2 (note that these values correspond also to those obtained with a mean value of the values for e = 10 -1 , e = 10 -2 and e = 10 -3 ). We observe that for N = 51 time windows, we gain approximately a factor 26.29 for Parareal and a factor 57.38 for Parareal-OSWRG which is better than the expected parallel efficiency. Case 2. We consider nonconforming meshes, refined in the region of the boudary layer (see Figure 6 (middle)), with a mesh size equal to 0.0065, 0.0131, and 0.0295 for the subdomains with a boundary along {x = 0}, those with boundaries along {x = 1 3 } and {x = 2 3 }, and those with a boundary along {x = 1}, respectively. The time step of the fine solver is ∆t F = 0.01. In Figure 7 (right), we plot the evolution of the relative error in L ∞ (0, T ; X )-norm, between the Parareal-OSWRG solution and the converged Parareal solution, as a function of the number of Parareal iterations, for different values of L. The Parareal algorithm (L = ∞) corresponds to the case L = 12. As for case 1, the horizontal dashed green lines represent three possible values of 10% of the error, and Table 10 shows the number of Parareal and total OSWRG iterations to reach these different values, for L = 1, 2, 4, 8, 12. We observe that L = 2 is the case with the fewest iterations. We also observe on this figure two convergence regimes as for Case 1.

The OSWRG algorithm used as a solver for problem (8.1) will need 24, 40 and 56 iterations to get a relative error (between the iterative solution and the converged OSWRG solution) smaller than 10 -2 , 10 -3 and 10 -4 , respectively. Table 11 shows the gain, when the fine solvers are performed in parallel, of the Parareal or the coupled Parareal-OSWRG methods compared to the OSWRG algorithm. In this table we take the values obtained for e = 10 -2 (note that these values correspond also to those obtained with a mean value of the values for e = 10 -1 , e = 10 -2 and e = 10 -3 ). We observe that for N = 51 time windows, the results with a conforming mesh (Table 9) or with nonconforming meshes (Table 11) are very close and slightly better for nonconforming meshes adapted to the physics. Indeed, the gain is approximately a factor 24.28 for Parareal and a factor 60 for Parareal-OSWRG which is significantly better than the expected parallel efficiency. In Figure 8, we plot the evolution of the relative error in L ∞ (0, T ; X )-norm, between the Parareal-OSWRG solution and the converged Parareal solution, as a function of the number of Parareal iterations, for different values of L. The Parareal algorithm (L = ∞) corresponds to the case L = 35. As for case 1, the horizontal dashed green lines represent three possible values of 10% of the error, and Table 12 shows the number of Parareal and total OSWRG iterations to reach 10 -3 for the values L = 1, 2, 4, 8, 35. As in Section 7.2, we observe that the pure Parareal algorithm doesn't perform well since it needs full 10 iterations to converge. Parareal-OSWRG with L = 2 is the best choice, and a good alternative method to the Parareal method. Indeed, as shown on Table 13, we gain a factor 29.58 for Parareal-OSWRG compared to the pure OSWRG algorithm, and a factor 14.58 compared to the pure Parareal algorithm. Remark 8.1. Note that in the results above, the number of OSWRG iterations to reach full convergence, inside each time window in the Parareal iterations [START_REF]An optimized Schwarz waveform relaxation method for the unsteady convection diffusion equation in two dimensions[END_REF], is smaller than the one needed to reach full convergence on (0, T ) (71), probably because the time windows are smaller. This explains the gain of a factor 2.03 when using Parareal, compared to the pure OSWRG algorithm.

Solver

Remark 8.2. The 2D simulations of this section don't show a comparison with the PSWR algorithm of [START_REF] Gander | A superlinear convergence estimate for the Parareal Schwarz waveform relaxation algorithm[END_REF]. However we have observed that the case L = 1 is not the optimal choice of L. Doing incomplete iterations for the coarse solver as well, will very likely lead to similar or slower convergence.

8.3.

Example in an industrial context. We consider a model problem given by ANDRA, the French National Agency for Radioactive Waste Management (see [START_REF] Hoang | Space-time domain decomposition methods for diffusion problems in mixed formulations[END_REF][START_REF] Omnes | A Posteriori Error Estimates and Stopping Criteria for Solvers using Domain Decomposition Method and with Local Time Stepping[END_REF][START_REF] Ali Hassan | A posteriori stopping criteria for space-time domain decomposition for the heat equation in mixed formulations[END_REF]), that simulates the transport of a contaminant in and around a nuclear waste repository site. The simulation domain is depicted in Figure 9 (left) (not at scale). The nuclear waste is stored in the repository (yellow), which is a 2950m by 10m rectangle located in the center of a clay domain of 3950m by 140m. In this example, we are concerned with the following time-dependent diffusion problem in mixed formulation:

σ σ σ = -S S S∇u in Ω × (0, T ), (8.3a)
φ ∂u ∂t + ∇ • σ σ σ = f in Ω × (0, T ), (8.3b)
where Ω = [0, 3950]×[0, 140] is decomposed into nine subdomains with Ω 5 the nuclear waste repository; see Figure 9. Here u is the (dimensionless) concentration of the contaminant, φ is the porosity, and S S S is the time-independent diffusion tensor, given in Table 14 (where I is the identity matrix). The initial condition is u 0 = 0, and we set homogeneous Dirichlet conditions on the top and the bottom of Ω, and homogeneous Neumann conditions on the other sides of ∂Ω. We are interested in the long-term behavior of the repository, over one million years, so that we set T = 10 6 years. A dimensionless form of this problem is given in [START_REF] Ali Hassan | A posteriori stopping criteria for space-time domain decomposition for the heat equation in mixed formulations[END_REF]. The number of triangles in the mesh of Ω is 58710 (see Figure 9 (right)) and the number of DoF is 148134. We use the Parareal-OSWRG method with Robin transmission conditions to solve this problem.

In Figure 10, we plot the relative error in L 2 (Ω × (0, T ))-norm, between the Parareal-OSWRG solution and the converged Parareal solution, versus the number of Parareal iterations, for different values of L. The Parareal algorithm (L = ∞) corresponds to L = 10.

In [START_REF] Omnes | A Posteriori Error Estimates and Stopping Criteria for Solvers using Domain Decomposition Method and with Local Time Stepping[END_REF], an a posteriori stopping criterion is given for this test case and the dashed green line shown in Figure 10 corresponds to this criterion. This line is reached after 3 iterations for Parareal and Parareal-OSWRG with L = 4, and after 13 iterations for Parareal-OSWRG with L = 2. 12 These data are reported in Table 15 (for the methods that converge in at most 30 iterations) and we observe that the fastest case is L = 4 with a total of 12 OSWRG iterations.

Table 16 shows the gain when the fine solvers are performed in parallel, of the Parareal or the coupled Parareal-OSWRG methods compared to the OSWRG method. In the first column we give the number of OSWRG iterations given in [START_REF] Omnes | A Posteriori Error Estimates and Stopping Criteria for Solvers using Domain Decomposition Method and with Local Time Stepping[END_REF] from the a posteriori stopping criterion, when the OSWRG method is used for solving problem (8.3). We see that for N = 10 time windows, we gain a factor 3.06 for Parareal and a factor 9.17 for Parareal-OSWRG which is almost the expected parallel efficiency. Remark 8.4. Note that to determine a suitable L in general, which may vary through Parareal iterations, one could design a posteriori error estimates that combine contributions from the scheme error and from the current nonconverged OSWR iteration. This would give a suitable L and an efficient stopping criterion for the Parareal-OSWR iterations, which would enhance the application of the method. 9. Proof of estimate (4.13) in Theorem 4.2 .

Proof. For i = 1, 2, let (u i , ξ i ) be defined by (4.4)-(4.5) with initial condition ū0 . We set ζ i := ξ i -ξ i , where ξ i is defined in (2.4), and e i := u i -u| Ωi , where u is the solution of (2.1). With (ζ 0 1 , ζ 0 2 ) given, the error e i satisfies for = 1, 2, . . . , L :

Le i = 0 in Ω i × I, (9.1a
(ζ i ) 2 = 1 2p i (ζ -1 i ) 2 .
Then, integrating on I, and considering L 2 (0, T ; H 1 (Ω i )) endowed with the norm u 2 L 2 (0,T ;H 1 (Ωi)) := I ν ∂ x u 2 + b u 2 dt, we obtain Summing over from 1 to L, we get a telescopic sum on the interface and thus we get (4.13).

Algorithm 3 . 1

 31 (Parareal) Choose an initial data (U 0 n ) n∈ 0,N with U 0 0 = u 0 and U 0 n an approximation of u(•, T n ), for example : U 0 n := G(I n-1 , U 0 n-1 ), for n = 1, 2, . . . , N . for k = 0, 1, . . . (Parareal iterations) do Set U k+1 0 = u 0 and perform the correction iterations

  ρ c (p) := max π T ≤ω≤ π ∆t |ρ 0 (p, ω)|, then the optimized Robin parameter p c is chosen such that it verifies ρ c (p c ) = min p>0 ρ c (p).(4.24)

. 3 )

 3 Update the interface term: ζ ζ ζ k+1,0 n = ζ ζ ζ k,L n . (6.4) end for The proof of Theorem 2.1 relies on the OSWR estimate (4.13) with e 0 = E k n (from incomplete iterations) and ζ ζ ζ = ζ ζ ζ k, n . Then summing with respect to k, the terms with ζ ζ ζ k, n are treated using the update ζ ζ ζ k+1,0 n = ζ ζ ζ k,L n

Fig. 1 .

 1 Fig. 1. Relative error versus Parareal iterations

  33 for Parareal-OSWR or PSWR.

Fig. 2 . 3

 23 Fig. 2. Relative error versus Parareal iterations

Fig. 3 .

 3 Fig. 3. Example 1: Rotating velocity (left), decomposition of Ω into four subdomains (middle), and decomposition of Ω into nine subdomains (right)

Fig. 4 . 1 : 6

 416 Fig. 4. Example 1: Relative error versus Parareal iterations. Case ν = 0.05 with four subdomains (left) and case ν = 0.1 with nine subdomains (right)

Fig. 5 .

 5 Fig. 5. Example 1 (case 1): Relative error versus Parareal iterations. Parareal-OSWRG with Ventcel (top left), Parareal-OSWR with Ventcel (top right), Parareal-OSWRG with Robin (bottom left), and Parareal-OSWR with Robin (bottom right)

  (left). The diffusion coefficient is ν = 0.01 and the final time is T = 51. The number of time windows for the Parareal iterations is N = 51. The time step of the coarse solver is ∆t C = 1.

Fig. 6 .

 6 Fig. 6. Example 2: Velocity field (left), nonconforming meshes (middle), and computed Parareal-OSWRG solution at final time (right)

Fig. 7 .

 7 Fig. 7. Example 2: Relative error versus Parareal iterations, with a zoom on the first iterations. Left : Case 1 (conforming meshes). Right : Case 2 (nonconforming meshes)

8

 8 Example 2 (case 1): Number of Parareal iterations k and total number of OSWRG iterations L * k, versus L to reach 3 different values: 10 -2 (top left), 10 -3 (top right) and 10 -4 (bottom)

L

  2): Gain factor for Parareal and coupled Parareal-OSWRG methods compared to the OSWRG solver for solving problem (8.1), with N = 51 time windows Case 3. We consider the example of case 1, now with T = 1 and N = 10, i.e. an advection dominated case on short time interval for which the Parareal algorithm converges slowly.

Fig. 8 .

 8 Fig. 8. Relative error versus Parareal iterations

  case 3): Gain factor for Parareal and coupled Parareal-OSWRG methods compared to the OSWRG solver for solving problem (8.1), with N = 10 time windows, to reach 10 -3

Fig. 9 . 2 ×

 92 Fig.9. Geometry of the nuclear waste repository (yellow) and the clay layer around it, with its decomposition into nine subdomains (left), and mesh used in and around the repository (right)

Fig. 10 .

 10 Fig. 10. Relative error versus Parareal iterations

Table 16 Example 3 :

 163 Gain factor for Parareal and coupled Parareal-OSWRG methods compared to the OSWRG solver for solving problem (8.3), with N = 10 time windows Remark 8.3. Note that increasing the contrast in the diffusion tensor S S S by taking S S S = 2 × 10 -9 in Ω 5 and S S S = 5 × 10 -13 in the other subdomains leads to slightly better results : The best choice is still L = 4, and for N = 10 time windows we gain a factor 2.67 for Parareal, and a factor 10 for Parareal-OSWRG, which is the expected parallel efficiency.

1 i 2 d dt e 1 2 + ν ∂ x e 1 2 + b e 1 2 -ν∂ x e 1 - a 2 e 1 2 d dt e 2 2 + ν ∂ x e 2 2 + b e 2 2 + 2 1 2 d dt e i 2 + ν ∂ x e i 2 + 1 2 d dt e i 2 + ν ∂ x e i 2 + b e i 2

 1222212222222212222 ) e i (•, 0) = e 0,i in Ω i , (9.1b)(B i e i )(0, •) = ζ -on I,(9.1c)whereζ i = B i e j (0, •), ≥ 1. (9.2)Then the proof follows the same steps as in[START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF] Theorem 5.15]. Multiplying equation(9.1a) by e i and integrating in space on Ω i , we obtain 1 (0, .)e 1 (0, .) = 0,1 ν∂ x e 2 -a2e 2 (0, .)e 2 (0, .) = 0.Introducing the Robin interface operators B i , i = 1, 2, defined in (2.2), and rewriting the terms on the interface in the formν∂ x e i -a 2 e i e i = 1 2p (B 1 e i ) 2 -(B 2 e i ) 2 , i = 1, 2,we obtain the energy estimates fori = 1, b e i 2 + 1 2p ((B j e i )(0, .)) 2 = 1 2p ((B i e i )(0, .)) 2 , j = 3 -i. Replacing (B i e i )(0, .) = ζ -1 i, (B j e i )(0, .) = ζ j (from (9.1c) and (9.2)), and summing the above expression over i, we get for all t ∈ I, i

1 2 ee i 2 L 2

 222 (., T ) 2 + i (0,T,H 1 (Ωi)) + 1 2p ζ ζ ζ 2 I = 1 2 i e i (., 0) 2 + 1 2p ζ ζ ζ -1 2 I .

  and there exists a constant C s.t.4

Table 1

 1 Number of Parareal iterations k and total number of OSWR iterations L * k, versus L

	L	PSWR 1	2	4	8	48 (∞)
	k	76	77 39 19 13	7
	L * k	77	77 78 76 104	336

Table 4

 4 Example in 1D (advection dominated): Gain factor for Parareal, coupled Parareal-OSWR, and PSWR methods compared to the OSWR solver, with N = 10 time windows

		OSWR Parareal Parareal-OSWR (L = 1) PSWR
	Iterations (n)	6	80	11	11
	Loss factor ( = n/6)	-	13.3	1.83	1.83
	Final gain factor (N/ )	-	0.75	5.45	5.45

Table 7

 7 Example 1 (case 2): Gain factor for Parareal and coupled Parareal-OSWRG methods compared to the OSWRG solver for solving problem (8.1), with N = 21 time windows

Table 9

 9 Example 2 (case 1): Gain factor for Parareal and coupled Parareal-OSWRG methods compared to the OSWRG solver for solving problem (8.1), with N = 51 time windows

Table 10

 10 Example 2 (case 2): Number of Parareal iterations k and total number of OSWRG iterations L * k, versus L to reach 3 different values: 10 -2 (top left), 10 -3 (top right) and 10 -4 (bottom)

		1	2	4	8 12 (≈ ∞)	L	1	2	4	8 12 (≈ ∞)
	k	16 6	4	3	3	k		47 17 11 7
	L * k 16 12 16 24	36	L * k 47 34 44 56	84
				L	1	2	4	8 12 (≈ ∞)
				k	79 28 17 11
				L * k 79 56 68 88	132
		Solver			OSWRG Parareal Parareal-OSWRG (L = 2)
		Iterations (n)			40	84		34
		Loss factor ( = n/40)	-	2.1		0.85
		Final gain factor (N/ )	-	24.28		60
						Table 11
		Example 2 (case

Table 12

 12 Number of Parareal iterations k and total number of OSWRG iterations L * k, versus L

Table 14

 14 Porosity φ, diffusion tensor S S S, and source term f

Table 15

 15 Number of Parareal iterations k and total number of OSWRG iterations L * k, versus L

See[START_REF]Non-homogeneous Boundary Value Problems and Applications[END_REF] Theorem 6.2 (extended to R in space), with m = 1 and r = 0.

See[START_REF]Non-homogeneous Boundary Value Problems and Applications[END_REF] Theorem 6.2 (extended to R in space), with m = m j = 1 and r = 0.

See[START_REF]Non-homogeneous Boundary Value Problems and Applications[END_REF] Theorem 2.1, with r = 2 and s = 1.

Note that we have a better result : u(•, T ) ∈ H 1 (Ω).

In the 1D case considered here, the interface is reduced to one point in space, thus the exchanged data depend on time only.

See[START_REF] Lions | Non-homogeneous Boundary Value Problems and Applications[END_REF] Theorem 11.4, with s = 1 4 .

Note that (4.18) is in the sense of distributions. Then (4.20) implies that (4.18) also holds in the classical sense.

Note that here the term "∂x" in B i is multiplied by ν while this is not the case in[START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF]. Thus ê i is slightly different here from the one of[START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF].

In the sense that L will be smaller than the number of iterations required for convergence.

This assumption is for example satisfied for the implicit Euler scheme, and is deduced from the regularity of the elliptic problem: U ∆T -ν∂xxU + a∂xU + bU = U 0 ∆T .

D.Q. BUI, C. JAPHET, Y. MADAY, AND P. OMNES

These optimized parameters are variable and calculated on each edge of the interfaces.

Acknowledgment. The authors warmly thank the referees for their comments that helped us to improve the results of this paper.

* This work was supported by the ANR project CINE-PARA under grant ANR-15-CE23-0019.