Laurence Halpern
email: halpern@math.univ-paris13.fr

Caroline Japhet
email: japhet@cscamm.umd.edu

AND J Ér

Émie Szeftel
email: jeremie.szeftel@ens.fr.

OPTIMIZED SCHWARZ WAVEFORM RELAXATION AND DISCONTINUOUS GALERKIN TIME STEPPING FOR HETEROGENEOUS PROBLEMS

Keywords: Coupling heterogeneous problems, domain decomposition, optimized Schwarz waveform relaxation, time discontinuous Galerkin, nonconforming grids, error analysis AMS subject classifications. 65 M 15, 65M50, 65M55

We design and analyze Schwarz waveform relaxation algorithms for domain decomposition of advection-diffusion-reaction problems with strong heterogeneities. These algorithms rely on optimized Robin or Ventcell transmission conditions, and can be used with curved interfaces. We analyze the semi-discretization in time with discontinuous Galerkin as well. We also show twodimensional numerical results using generalized mortar finite elements in space.

1. Introduction. In many fields of applications such as reactive transport, far field simulations of underground nuclear waste disposal or ocean-atmosphere coupling, models have to be coupled in different spatial zones, with very different space and time scales and possible complex geometries. For such problems with long time computations, a splitting of the time interval into windows is essential, with the possibility to use robust and fast solvers in each time window.

The Optimized Schwarz Waveform Relaxation (OSWR) method was introduced for linear parabolic and hyperbolic problems with constant coefficients in [START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF]. It was analyzed for advection diffusion equations, and applied to non constant advection, in [START_REF] Martin | An optimized Schwarz waveform relaxation method for unsteady convection diffusion equation[END_REF]. The algorithm computes independently in each subdomain over the whole time interval, exchanging space-time boundary data through optimized transmission operators. The operators are of Robin or Ventcell type, with coefficients optimizing a convergence factor, extending the strategy developed by F. Nataf and coauthors [START_REF] Charton | Méthode de décomposition de domaine pour l'équation d'advection-diffusion[END_REF][START_REF] Japhet | Méthode de décomposition de domaine et conditions aux limites artificielles en mécanique des fluides : méthode Optimisée d'Ordre 2 (OO2)[END_REF]. The optimization problem was analyzed in [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF], [START_REF] Bennequin | A homographic best approximation problem with application to optimized Schwarz waveform relaxation[END_REF].

This method potentially applies to different space-time discretization in subdomains, possibly nonconforming and needs a very small number of iterations to converge. Numerical evidences of the performance of the method with variable smooth coefficients were given in [START_REF] Martin | An optimized Schwarz waveform relaxation method for unsteady convection diffusion equation[END_REF]. An extension to discontinuous coefficients was introduced in [START_REF] Gander | Schwarz waveform relaxation method for advectiondiffusion-reaction problems with discontinuous coefficients and non-matching grids[END_REF], with asymptotically optimized Robin transmission conditions in some particular cases.

The discontinuous Galerkin finite element method in time offers many advantages. Rigorous analysis can be made for any degree of accuracy and local time-stepping, and time steps can be adaptively controlled by a posteriori error analysis, see [START_REF] Thomee | Galerkin Finite Element Methods for Parabolic Problems[END_REF][START_REF] Johnson | Time discretization of parabolic problems by the discontinuous Galerkin method[END_REF][START_REF] Makridakis | A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF]. In a series of presentations in the regular domain decomposition meeting we presented the DG-OSWR method, using discontinuous Galerkin for the time discretization of the OSWR. In [START_REF] Blayo | Optimized Schwarz waveform relaxation algorithms with nonconforming time discretization for coupling convection-diffusion problems with discontinuous coefficients[END_REF], [START_REF] Halpern | Discontinuous Galerkin and Nonconforming in Time Optimized Schwarz Waveform Relaxation for Heterogeneous Problems[END_REF], we introduced the algorithm in one dimension with discontinuous coefficients. In [START_REF] Halpern | Discontinuous Galerkin and nonconforming in time optimized Schwarz waveform relaxation[END_REF], we extended the method to the two dimensional case. For the space discretization, we extended numerically the nonconforming approach in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF] to advection-diffusion problems and optimized order 2 transmission conditions, to allow for non-matching grids in time and space on the boundary. The space-time projections between subdomains were computed with an optimal projection algorithm without any additional grid, as in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF]. Two dimensional simulations were presented. In [START_REF] Halpern | Space-time nonconforming Optimized Schwarz Waveform Relaxation for heterogeneous problems and general geometries[END_REF] we extended the proof of convergence of the OSWR algorithm to nonoverlapping subdomains with curved interfaces. Only sketches of proofs were presented.

The present paper intends to give a full and self-contained account of the method for the advection diffusion reaction equation with non constant coefficients.

In Section 2, we present the Robin and Ventcell algorithms at the continuous level in any dimension, and give in details the new proofs of convergence of the algorithms for nonoverlapping subdomains with curved interfaces. Then in Section 3, we discretize in time with discontinuous Galerkin, and prove the convergence of the semi-discrete algorithms for flat interfaces. Error estimates are derived from the classical ones [START_REF] Thomee | Galerkin Finite Element Methods for Parabolic Problems[END_REF]. The fully discrete problem is introduced in Section 4, using finite elements. The interfaces are treated by a new cement approach, extending the method in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF]. Given the length of the paper, the numerical analysis will be treated in a forthcoming paper. We finally present in Section 5 simulations for two subdomains, with piecewise smooth coefficients and a curved interface, for which no error estimates are available. We also include an application to the porous media equation.

Consider the advection-diffusion-reaction equation in

Ω = R N ∂ t u + ∇ • (b b bu -ν∇u) + cu = f in Ω × (0, T ), (1.1)
with initial condition u(0, x) = u 0 (x) x ∈ Ω.

(1.

2)

The advection and diffusion coefficients b b b and ν , as well as the reaction coefficient c, are piecewise smooth, the problem is parabolic, i.e. ν ≥ ν 0 > 0 a.e. in R N . Theorem 1.1 (Well-posedness and regularity, [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]). Let Ω = R N . Suppose b b b ∈ (W 1,∞ (Ω)) N , ν ∈ W 1,∞ (Ω) and c ∈ L ∞ (Ω). If the initial value u 0 is in H 1 (Ω), and the right-hand side f is in L 2 (0, T ; L 2 (Ω)), then there exists a unique solution u of (1.1), (1.2) in H 1 (0, T ; L 2 (Ω)) ∩ L ∞ (0, T ; H 1 (Ω)) ∩ L 2 (0, T ; H 2 (Ω)).

We consider now a decomposition of Ω into nonoverlapping subdomains Ω i , i ∈ [ [1, I]], as depicted in Figure 1.1. In all cases the boundaries between the subdomains are supposed to be hyperplanes at infinity. Problem (1.1) is equivalent to solving I problems in subdomains Ω i , with transmission conditions on the interface Γ i,j between two neighboring subdomains Ω i and Ω j , given by the jumps [u] = 0, [(ν∇ub b b) • n n n i ] = 0. Here n n n i is the unit exterior normal to Ω i . As coefficients ν and b b b are possibly discontinuous on the interface, we note, for s ∈ Γ i,j , ν i (s) = lim ε→0+ ν(sεn n n i ). The same notation holds for b b b and c.

Ω i Ω i
To any i ∈ [ [1, m]], we associate the set N i of indices of the neighbors of Ω i . Following [START_REF] Charton | Méthode de décomposition de domaine pour l'équation d'advection-diffusion[END_REF][START_REF] Gander | Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation[END_REF][START_REF] Gander | Schwarz waveform relaxation method for advectiondiffusion-reaction problems with discontinuous coefficients and non-matching grids[END_REF], we propose as preconditioner for (1.1, 1.2), the sequence of problems

∂ t u k i + ∇ • (b b b i u k i -ν i ∇u k i ) + c i u k i = f in Ω i × (0, T ), (1.3a
)

ν i ∂ n n n i -b b b i • n n n i u k i + S i,j u k i = ν j ∂ n n n i -b b b j • n n n i u k-1 j + S i,j u k-1 j on Γ i,j , j ∈ N i . (1.3b)
The boundary operators S i,j , acting on the part Γ i,j of the boundary of Ω i shared by the boundary of Ω j are given by

S i,j ϕ = p i,j ϕ + q i,j (∂ t ϕ + ∇ Γi,j • (r r r i,j ϕ -s i,j ∇ Γi,j ϕ)). (1.4)
∇ Γ and ∇ Γ • are respectively the gradient and divergence operators on Γ. p i,j , q i,j , s i,j are functions in L ∞ (Γ i,j ) and r r r i,j is in (L ∞ (Γ i,j )) N -1 . The initial value is that of u 0 in each subdomain. An initial guess (g i,j ) is given on

L 2 ((0, T ) × Γ i,j ) for i ∈ [[1, I]], j ∈ N i .
By convention for the first iterate, the right-hand side in (1.3) is given by g i,j . Under regularity assumptions, solving (1.1) is equivalent to solving

∂ t u i + ∇ • (b b b i u i -ν i ∇u i ) + c i u i = f in Ω i × (0, T ), ν i ∂ n n n i -b b b i • n n n i u i + S i,j u i = ν j ∂ n n n i -b b b j • n n n i u j + S i,j u j on Γ i,j × (0, T ), j ∈ N i , (1.5
) for i ∈ [[1, I]] with u i the restriction of u to Ω i .
2. Studying the algorithm for the P.D.E. The first step of the study is to give a frame for the definition of the iterates.

∂ t w + ∇ • (b b bw -ν∇w) + cw = f in O × (0, T ), ν ∂ n n n w -b b b • n n n w + Sw = g on Γ × (0, T ), w(•, 0) = u 0 in O, (2.1)
where n n n is the exterior unit normal to O. The boundary operator S is defined on Γ = ∂O by

Sw = p w + q (∂ t w + ∇ Γ • (r r rw -s∇ Γ w)).
(2.

2)

The domain O has either form depicted in Figure 1.1, left for q = 0, right otherwise. The functions p, q and s are in L ∞ (Γ), and r r r is in (L ∞ (Γ)) N -1 . Either q = 0, and the boundary condition is of Robin type, or we suppose q ≥ q 0 > 0 and the operator will be referred to as Order 2 or Ventcell operator. In the latter case, we need the spaces

H s s (O) = {v ∈ H s (O), v |Γ ∈ H s (Γ)}, (2.3)
which are defined for s > 1/2, and equipped with the scalar product

(w, v) H s s (O) = (w, v) H s (O) + (qw, v) H s (Γ) . (2.4)
We define the bilinear forms m on H 1 (O) and a on

H 1 1 (O) by m(w, v) = (w, v) L 2 (O) + (qw, v) L 2 (Γ) , (2.5)
and

a(w, v) := O ( 1 2 ((b b b • ∇w)v -(b b b • ∇v)w)) dx + O ν∇w • ∇v dx + O (c + 1 2 ∇ • b b b)wv dx + Γ (p - b b b • n n n 2 + q 2 ∇ Γ • r r r)wv + q 2 (∇ Γ • (r r rw)v -∇ Γ • (r r rv)w) + qs∇ Γ w • ∇ Γ v dσ, (2.6)
By the Green's formula, we can write a variational formulation of (2.1):

d dt m(w, v) + a(w, v) = (f, v) L 2 (O) + (g, v) L 2 (Γ) . (2.7)
The well-posedness is a generalization of results in [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Bennequin | A homographic best approximation problem with application to optimized Schwarz waveform relaxation[END_REF][START_REF] Szeftel | Calcul pseudo-différentiel et para-différentiel pour l'étude des conditions aux limites absorbantes et des propriétés qualitatives des EDP non linéaires[END_REF]. It relies on energy estimates and Grönwall's lemma.

Theorem 2.1. Suppose ν ∈ W 1,∞ (O), b b b ∈ (W 1,∞ (O)) N , c ∈ L ∞ (O), p ∈ L ∞ (Γ), q ∈ L ∞ (Γ), r r r ∈ (W 1,∞ (Γ)) N -1 , s ∈ W 1,∞ (Γ) with s > 0 a.e. If q = 0, if f is in H 1 (0, T ; L 2 (O)), u 0 is in H 2 (O) and g is in H 1 (0, T ; L 2 (Γ)) ∩ L ∞ (0, T ; H 1/2 (Γ)), satisfying the compatibility condition ν ∂ n n n u 0 -b b b • n n n u 0 + pu 0 = g , the subdomain problem (2.1) has a unique solution w in L ∞ (0, T ; H 2 (O)) ∩ W 1,∞ (0, T ; L 2 (O)). If q ≥ q 0 > 0 a.e., if f is in H 1 (0, T ; L 2 (O)), u 0 is in H 2 2 (O), and g is in H 1 ((0, T ); L 2 (Γ)), problem (2.1) has a unique solution w in L ∞ (0, T ; H 2 2 (O)) ∩ W 1,∞ (0, T ; L 2 (O)) with ∂ t w ∈ L ∞ (0, T ; L 2 (Γ)).
Proof. The existence result relies on a Galerkin method like in [START_REF] Szeftel | Absorbing boundary conditions for reaction-diffusion equations[END_REF][START_REF] Szeftel | Calcul pseudo-différentiel et para-différentiel pour l'étude des conditions aux limites absorbantes et des propriétés qualitatives des EDP non linéaires[END_REF]. In the sequel, α, β, • • • denote positive real numbers depending only on the coefficients and the geometry. The basic estimate is obtained by multiplying the equation by w and integrating by parts in the domain. We set w = w L 2 (O) and w Γ = w L 2 (Γ) . 1 2

d dt m(w, w) + a(w, w) = (f, w) L 2 (O) + (g, w) L 2 (Γ) .
With the assumptions on the coefficients, we have

Case q = 0. a(w, w) = O ν|∇w| 2 dx + O (c + 1 2 ∇ • b b b)w 2 dx + Γ (p - b b b • n n n 2 w 2 dσ ≥ α( ∇w 2 -β( w 2 + w 2 Γ ) ≥ α 2 ∇w 2 -γ w 2 ,
the last inequality coming from the trace theorem

w 2 Γ ≤ C ∇w w . (2.8)
We obtain with the Cauchy-Schwarz inequality

1 2 d dt w 2 + α 4 ∇w 2 ≤ η w 2 + δ( f 2 + g 2 Γ )
.

We now have with Grönwall's lemma

w(t) 2 + α 2 t 0 ( ∇w(s) 2 ds ≤ e 2ηT ( u 0 2 + 2δ( f 2 L 2 (0,T ;L 2 (O)) + g 2 L 2 (0,T ;L 2 (Γ)) )).
(2.9)

We apply (2.9) to w t :

w t (t) 2 + α 2 t 0 ∇w t (s) 2 ds ≤ e 2ηT ( w t0 2 + 2δ( f t 2 L 2 (0,T ;L 2 (O)) + g t 2 L 2 (0,T ;L 2 (Γ)) )).
Thanks to the compatibility condition, w t0 can be estimated, using the equation, by

w t0 ≤ ζ( u 0 H 2 (O) + f (0, •)
) , and we obtain

w t (t) 2 + α 2 t 0 ∇w t (s) 2 ds ≤ θe 2ηT ( u 0 H 2 (O) + ( f 2 H 1 (0,T ;L 2 (O)) + g 2 H 1 (0,T ;L 2 (Γ)) )).
Case q ≥ q 0 > 0 a.e.

a(w, w) = O ν|∇w| 2 dx + O (c + 1 2 ∇ • b b b)w 2 dx + Γ (p - b b b • n n n 2 + q 2 ∇ Γ • r r r)w 2 + qs|∇ Γ w| 2 dσ, ≥ α( ∇w 2 + ∇ Γ w 2 Γ ) -βm(w, w),
and by the Cauchy-Schwarz inequality

O (f, w) dx + Γ (g, w) Γ dσ ≤ γm(w, w) + δ( f 2 + g 2 Γ ).
Collecting these inequalities, we obtain

1 2 d dt m(w, w) + α( ∇w 2 + ∇ Γ w 2 Γ ) ≤ (β + γ)m(w, w) + δ( f 2 + g 2 Γ ).
We now integrate in time and use Grönwall's lemma to obtain for any t in (0, T )

m(w(t), w(t)) + 2α t 0 ( ∇w(s) 2 + ∇ Γ w(s) 2 Γ )ds ≤ e (β+γ)T (m(u 0 , u 0 ) + 2δ( f 2 L 2 (0,T ;L 2 (O)) + g 2 L 2 (0,T ;L 2 (Γ)) )).
(2.10)

We apply (2.10) to w t :

m(w t (t), w t (t)) + 2α t 0 ( ∇w t (s) 2 + ∇ Γ w t (s) 2 Γ )ds ≤ e (β+γ)T (m(w t0 , w t0 ) + 2δ( f t 2 L 2 (0,T ;L 2 (O)) + g t 2 L 2 (0,T ;L 2 (Γ)) )).
We now use the equations at time 0 to estimate m(w t0 , w t0 ). From the equation in the domain, we deduce that

w t0 ≤ ζ( u 0 H 2 (O) + f (0, •) ),
and from the boundary condition that

w t0 Γ ≤ η( u 0 H 2 2 (O) + g(0, •) Γ ), which gives altogether m(w t (t), w t (t)) + 2α t 0 ( ∇w t (s) 2 + ∇ Γ w t (s) 2 Γ )ds ≤ θe (β+γ)T ( u 0 2 H 2 2 (O) + f 2 H 1 (0,T ;L 2 (O)) + g 2 H 1 (0,T ;L 2 (Γ)) )).
(2.11)

We can now apply the Galerkin method. When q = 0, we work in

H 1 (0, T ; H 1 (O)) ∩ W 1,∞ (0, T ; L 2 (O)) , while if q ≥ q 0 > 0 a.e we consider H 1 (0, T ; H 1 1 (O)) ∩ W 1,∞ (0, T ; L 2 (O)) ∩ W 1,∞ (0, T ; L 2 (Γ))
. This gives a unique solution w. The regularity H 2 is obtained for q = 0 by the usual regularity results for the Laplace equation with Neumann boundary condition, since

-∆w = 1 ν (f -w t -∇ • (b b bw) + ∇ν • ∇w) ∈ L ∞ (0, T ; L 2 (O)), ∂ n n n w = 1 ν (b b b • n n n -p)w + 1 ν g ∈ L ∞ (0, T ; H 1/2 (Γ)).
In the other case, we have that

-∆w = 1 ν (f -w t -∇ • (b b bw) + ∇ν • ∇w) ∈ L ∞ (0, T ; L 2 (O)), ν∂ n n n -qs∆ Γ w = 1 ν (b b b • n n n -p -q (∂ t + ∇ Γ • r r r -(∇ Γ • s)∇ Γ )))w + 1 ν g ∈ L ∞ (0, T ; L 2 (Γ)),
and we conclude like in [START_REF] Szeftel | Absorbing boundary conditions for reaction-diffusion equations[END_REF][START_REF] Szeftel | Calcul pseudo-différentiel et para-différentiel pour l'étude des conditions aux limites absorbantes et des propriétés qualitatives des EDP non linéaires[END_REF].

2.2. Convergence analysis for Robin transmission conditions. We suppose here the coefficients q i to be zero everywhere. Given initial guess (g i,j ) on L 2 ((0, T ) × Γ i,j ) for i ∈ [[1, I]], j ∈ N i , the algorithm reduces in each subdomain to

∂ t u k i + ∇ • (b b b i u k i -ν i ∇u k i ) + c i u k i = f in Ω i × (0, T ), (2.12a
)

ν i ∂ n n ni -b b b i • n n n i u k i + p i,j u k i = ν j ∂ n n ni -b b b j • n n n i u k-1 j + p i,j u k-1 j on Γ i,j , j ∈ N i . (2.12b)
The well-posedness for the boundary value problem in the previous section permits to define the sequence of iterates. We now consider the convergence of this sequence.

Theorem 2.2. For coefficients p i,j such that p i,j + p j,i > 0 a.e., the sequence (u k i ) k∈N of solutions of (2.12) converges to the solution u of problem (1.1). Proof. By linearity, it is sufficient to prove that the sequence of iterates converges to zero if f = u 0 = 0.

We multiply (2.12a) by u k i , integrate on Ω i , and use the Green's formula. We obtain 1 2

d dt u k i (t, •) 2 L 2 (Ωi) + (ν i ∇u k i , ∇u k i ) L 2 (Ωi) + ((c i + 1 2 ∇ • b b b i )u k i , u k i ) L 2 (Ωi) - j∈Ni Γi,j ν i ∂ n n n i u k i - b b b i • n n n i 2 u k i u k i dσ = 0. (2.13)
We use now

ν i ∂ n n ni u k i -b b b i • n n n i u k i + p i,j u k i 2 -ν i ∂ n n ni u k i -b b b i • n n n i u k i -p j,i u k i 2 = 2(p i,j + p j,i ) ν i ∂ n n n i u k i - b b b i • n n n i 2 u k i u k i + (p i,j + p j,i )(p i,j -p j,i -b b b i • n n n i )(u k i ) 2 . (2.14)
We replace the boundary term in (2.13), and integrate in time. Since the initial value vanishes, we have for any time t,

u k i (t) 2 L 2 (Ωi) + 2 t 0 (ν i ∇u k i , ∇u k i ) L 2 (Ωi) + ((c i + 1 2 ∇ • b b b i )u k i , u k i ) L 2 (Ωi) dτ + j∈Ni t 0 Γi,j 1 p i,j + p j,i ν i ∂ n n n i u k i -b b b i • n n n i u k i -p j,i u k i 2 dσ dτ = j∈Ni t 0 Γi,j 1 p i,j + p j,i ν i ∂ n n n i u k i -b b b i • n n n i u k i + p i,j u k i 2 dσ dτ + j∈Ni t 0 Γi,j (p j,i -p i,j -b b b i • n n n i )(u k i ) 2 dσ dτ.
Since the coefficients are all bounded, the last term in the right-hand side can be handled by the trace theorem (2.8) to be canceled with the terms in the left-hand side like in the proof of Theorem 2.2. We further insert the transmission condition in the right-hand side:

u k i (t) 2 L 2 (Ωi) +ν 0 t 0 ∇u k i 2 L 2 (Ωi) dτ + j∈Ni t 0 Γi,j 1 p i,j + p j,i ν i ∂ n n n i u k i -b b b i •n n n i u k i -p j,i u k i 2 dσ dτ ≤ j∈Ni t 0 Γi,j 1 p i,j + p j,i (ν j ∂ n n n i u k-1 j -b b b j •n n n i u k-1 j +p i,j u k-1 j 2 dσ dτ +C 1 t 0 u k i 2 L 2 (Ωi) .
We sum on the subdomains, and on the iterations, the boundary terms cancel out except the first and last ones, and we obtain for any t ∈ (0, T ),

k∈[[1,K]] i∈[[1,I]] u k i (t) 2 L 2 (Ωi) +ν 0 t 0 ∇u k i 2 L 2 (Ωi) dτ ≤ α(t)+C 1 k∈[[1,K]] i∈[[1,I]] t 0 u k i 2 L 2 (Ωi) , (2.15)
with

α(t) = i∈[[1,I]] j∈Ni t 0 Γi,j 1 p i,j + p j,i (ν j ∂ n n ni u 0 j -b b b j • n n n i u 0 j + p i,j u 0 j 2 dσ dτ.
We now apply Grönwall's lemma and obtain that for any K > 0,

k∈[[1,K]] i∈[[1,I]] u k i (t) 2 L 2 (Ωi) ≤ α(T )e C1T ,
which proves that the sequence u k i converges to zero in L 2 ((0, T ) × Ω i ) for each i, and concludes the proof of the theorem.

Remark 2.3. [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF]) and we conclude without using Grönwall's lemma.

In the case ∇ • b b b = 0, if p j,i -p i,j -b b b i • n n n i = 0 and c i ≥ α 0 > 0, then C 1 = 0 in (2.

Order 2 transmission conditions.

Theorem 2.4. Assume p i,j ∈ W 1,∞ (Ω i ), p i,j + p j,i > 0 a.e., q i,j

= q > 0, b b b i ∈ (W 1,∞ (Ω i )) N , ν i ∈ W 1,∞ (Ω), r r r i,j ∈ (W 1,∞ (Ω i )) N -1 ,
with r r r i,j = r r r j,i on Γ i,j , s i,j ∈ W 1,∞ (Ω i ), s i,j > 0 with s i,j = s j,i on Γ i,j , and the domain is cut in bands as in Figure 1.1, right. Then, the algorithm (1.3) converges in each subdomain to the solution u of problem (1.1).

Proof. We first need some results in differential geometry. For any i ∈ [[1, I]], For every j ∈ N i , the normal vector n n n i can be extended in a neighbourhood of Γ i,j in Ω i as a smooth function ñ n n i with length one. Let ψ i,j ∈ C ∞ (Ω i ), such that ψ i,j ≡ 1 in a neighbourhood of Γ i,j , ψ i,j ≡ 0 in a neighbourhood of Γ i,k for k ∈ N i , k = j and j∈Ni ψ i,j > 0 on Ω i . We can assume that ñ n n i is defined on a neighbourhood of the support of ψ i,j . We extend the tangential gradient and divergence operators in the support of ψ i,j by:

∇ Γi,j ϕ := ∇ϕ -(∂ ñ n ni ϕ)ñ n n i , ∇ Γi,j • ϕ ϕ ϕ := ∇ • (ϕ ϕ ϕ -(ϕ ϕ ϕ • ñ n n i )ñ n n i ).
It is easy to see that ( ∇ Γi,j ϕ) |Γi,j = ∇ Γi,j ϕ, ( ∇ Γi,j • ϕ ϕ ϕ) |Γi,j = ∇ Γi,j • ϕ ϕ ϕ and for ϕ ϕ ϕ and χ with support in supp(ψ i,j ), we have

Ωi ( ∇ Γi,j • ϕ ϕ ϕ) χ dx = - Ωi ϕ ϕ ϕ • ∇ Γi,j χ dx. (2.16)
Now we prove Theorem 2.4. We consider the algorithm (1.3) on the error, so we suppose f = u 0 = 0. We set

ϕ i = ϕ L 2 (Ωi) , ϕ 2 i = √ ν i ∇ϕ 2 i , ϕ i,∞ = ϕ L ∞ (Ωi) , ϕ i,1,∞ = ϕ W 1,∞ (Ωi) and β i,j = pi,j +pj,i 2 .
The proof is based on energy estimates containing the term

t 0 Γi,j ν i ∂ n n n i u k i -b b b i • n n n i u k i + S i,j u k i 2 dσ dτ,
and that we derive by multiplying successively the first equation of (1.3) by the terms

β 2 i u k i , ∂ t u k i , ∇ Γi,j • (ψ 2 i,j r r r i,j u k i ) and -∇ Γi,j • (ψ 2 i,j s i,j ∇ Γi,j u k i ).
We multiply the first equation of (1.3) by β 2 i u k i , integrate on (0, t) × Ω i and integrate by parts in space,

1 2 β i u k i (t) 2 i + t 0 β i u k i (τ, •) 2 i dτ - t 0 Ωi β i (b b b i • ∇β i )(u k i ) 2 dx dτ + t 0 Ωi (c i + 1 2 ∇ • b b b i )β 2 i (u k i ) 2 dx dτ - t 0 Ωi ν i |∇β i | 2 (u k i ) 2 dx dτ - t 0 Γi,j β 2 i,j (ν i ∂ n n n i u k i - b b b i • n n n i 2 u k i ) u k i dσ dτ = 0. (2.17)
We multiply the first equation of (1.3) by ∂ t u k i , integrate on (0, t) × Ω i and then integrate by parts in space,

t 0 ∂ t u k i 2 i dτ + 1 2 u k i (t) 2 i + t 0 Ωi (c i u k i +∇•(b b b i u k i )) ∂ t u k i dx dτ - t 0 Γi,j ν i ∂ n n n i u k i ∂ t u k i dσ dτ = 0. (2.18)
We multiply the first equation of (1.3) by ∇ Γi,j • (ψ 2 i,j r r r i,j u k i ) integrate on (0, t) × Ω i and integrate by parts in space:

t 0 Ωi ∂ t u k i ∇ Γi,j • (ψ 2 i,j r r r i,j u k i ) dx dτ + t 0 Ωi ∇ • (b b b i u k i ) ∇ Γi,j • (ψ 2 i,j r r r i,j u k i ) dx dτ + t 0 Ωi c i u k i ∇ Γi,j • (ψ 2 i,j r r r i,j u k i ) dx dτ + t 0 Ωi ν i ∇u k i •∇ ∇ Γi,j •(ψ 2 i,j r r r i,j u k i ) dx dτ - t 0 Γi,j ν i ∂ n n n i u k i ∇ Γi,j •(r r r i,j u k i ) dσ dτ = 0. (2.19)
We observe that

t 0 Ωi ν i ∇u k i • ∇ ∇ Γi,j • (ψ 2 i,j r r r i,j u k i ) dx dτ = t 0 Ωi ν i ∇u k i •∇( ∇ Γi,j •(ψ 2 i,j r r r i,j )u k i ) dx dτ + t 0 Ωi ν i ∇u k i •∇(ψ 2 i,j r r r i,j • ∇ Γi,j u k i ) dx dτ, (2.20)
with t 0 Ωi ν i ∇u k i • ∇(ψ 2 i,j r r r i,j • ∇ Γi,j u k i ) dx dτ ≥ - 1 4 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i 2 i dτ -C t 0 Ωi ( ∇u k i 2 i + u k i 2 i ) dx dτ. (2.
t 0 Ω i ∂tu k i ∇Γ i,j • (ψ 2 i,j r r ri,j u k i ) + ∇ • (b b biu k i ) ∇Γ i,j • (ψ 2 i,j r r ri,j u k i ) + ciu k i ∇Γ i,j • (ψ 2 i,j r r ri,j u k i ) dx dτ - 1 4 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i ) 2 i dτ - t 0 Γi,j ν i ∂ n n ni u k i ∇ Γi,j • (r r r i,j u k i ) dσ dτ ≤ C t 0 Ωi ( √ ν i ∇u k i 2 i + β i u k i 2 i ) dx dτ.
(2.22) Now we multiply the first equation of (1.3) by -∇ Γi,j • (ψ 2 i,j s i,j ∇ Γi,j u k i ) integrate on (0, t) × Ω i and integrate by parts in space:

1 2 ψ i,j √ s i,j ∇ Γi,j u k i (t) 2 i - t 0 Ωi ∇ • (b b b i u k i ) ∇ Γi,j • (ψ 2 i,j s i,j ∇ Γi,j u k i ) dx dτ + t 0 Ω i ψ 2 i,j si,j ∇Γ i,j (ciu k i ) • ∇Γ i,j u k i dx dτ - t 0 Ω i νi∇u k i • ∇( ∇Γ i,j • (ψ 2 i,j si,j ∇Γ i,j u k i )) dx dτ + t 0 Γi,j ν i ∂ n n n i u k i ∇ Γi,j • (s i,j ∇ Γi,j u k i ) dσ dτ = 0. (2.23)
We have,

- t 0 Ωi ν i ∇u k i • ∇( ∇ Γi,j • (ψ 2 i,j s i,j ∇ Γi,j u k i )) dx dτ ≥ 1 2 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i ) 2 i dτ -C 1 t 0 ∇u k i 2 i dτ. (2.24)
Replacing (2.24) in (2.23) leads to

1 2 ψ i,j √ s i,j ∇ Γi,j u k i (t) 2 i + 1 2 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i ) 2 i dτ + t 0 Ωi ψ 2 i,j s i,j c i | ∇ Γi,j u k i | 2 dx dτ + t 0 Γi,j ν i ∂ n n ni u k i ∇ Γi,j • (s i,j ∇ Γi,j u k i ) dσ dτ ≤ t 0 Ωi ∇ • (b b b i u k i ) ∇ Γi,j • (ψ 2 i,j s i,j ∇ Γi,j u k i ) dx dτ + C t 0 √ ν i ∇u k i 2 i dτ. (2.25)
Multiplying (2.18), (2.22) and (2.25) by q, and adding the three equations with (2.17), we get

1 2 β i u k i (t) 2 i + q u k i (t) 2 i +q ψ i,j √ s i,j ∇ Γi,j u k i (t) 2 i + t 0 β i u k i (τ, •) 2 i dτ + q t 0 ∂ t u k i 2 i dτ + q 4 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i 2 i dτ - t 0 Γi,j β 2 i (ν i ∂ n n n i u k i - b b b i • n n n i 2 u k i ) u k i dσ dτ -q t 0 Γi,j ν i ∂ n n n i u k i ∂ t u k i + ∇ Γi,j • (r r r i,j u k i ) -∇ Γi,j • (s i,j ∇ Γi,j u k i ) dσ dτ ≤ q( 1 2 b b b i i,1,∞ + r r r i,j i,1,∞ ) t 0 u k i i ∂ t u k i i dτ + q( b b b i i,∞ + r r r i,j i,∞ ) t 0 ∇u k i i ∂ t u k i i dτ + q 2 b b b i i,∞ t 0 ∇u k i i ∇ Γi,j • (ψ 2 i,j s i,j ∇ Γi,j u k i ) i dτ + q( 1 2 b b b i i,1,∞ + c i i,∞ ) t 0 u k i i ∇ Γi,j • (ψ 2 i,j s i,j ∇ Γi,j u k i ) i dτ + q 2 ( b b b i i,1,∞ + c i i,∞ ) u k i (t) 2 i + C t 0 β i u k i 2 i dτ + q t 0 √ ν i ∇u k i 2 i dτ .
We bound the right-hand side by 1 2

q 2 t 0 ∂ t u k i 2 i dτ + q 4 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i 2 i dτ + q 2 ( b b b i i,1,∞ + c i i,∞ ) u k i (t) 2 i + C t 0 β i u k i 2 i dτ + q t 0 √ ν i ∇u k i 2 i dτ .
We simplify the terms which appear on both sides, and obtain

1 2 β i u k i (t) 2 i + q u k i (t) 2 i +q ψ i,j √ s i,j ∇ Γi,j u k i (t) 2 i + t 0 β i u k i (τ, •) 2 i dτ + q 2 t 0 ∂ t u k i 2 i dτ + q 8 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i 2 i dτ - t 0 Γi,j β 2 i (ν i ∂ n n ni u k i - b b b i • n n n i 2 u k i ) u k i dσ dτ -q t 0 Γi,j ν i ∂ n n n i u k i ∂ t u k i + ∇ Γi,j • (r r r i,j u k i ) -∇ Γi,j • (s i,j ∇ Γi,j u k i ) dσ dτ ≤ C t 0 β i u k i 2 i dτ + q t 0 √ ν i ∇u k i 2 i dτ . (2.26)
Recalling that s i,j = s j,i on Γ i,j and r r r i,j = r r r j,i on Γ i,j , we use now:

ν i ∂ n n ni u k i -b b b i • n n n i u k i + S i,j u k i 2 -ν i ∂ n n ni u k i -b b b i • n n n i u k i -S j,i u k i 2 = 4 β 2 i,j (νi∂n n n i u k i - b b bi • n n ni 2 u k i )u k i + qνi∂n n n i u k i (∂tu k i + ∇Γ i,j • (r r ri,j u k i ) -∇Γ i,j • (si,j ∇Γ i,j u k i )) + 2q(p i,j -p j,i -2b b b i • n n n i )(∂ t u k i + ∇ Γi,j • (r r r i,j u k i ) -∇ Γi,j • (s i,j ∇ Γi,j u k i ))u k i + (p i,j + p j,i )(p i,j -p j,i -b b b i • n n n i )(u k i ) 2 . (2.27)
Replacing (2.27) into (2.26), we obtain

1 2 β i u k i (t) 2 i + q u k i (t) 2 i +q ψ i,j √ s i,j ∇ Γi,j u k i (t) 2 i + t 0 β i u k i (τ, •) 2 i dτ + q 2 t 0 ∂ t u k i 2 i dτ + 1 4 t 0 Γi,j ν i ∂ n n n i u k i -b b b i • n n n i u k i -S j,i u k i 2 dσ dτ + q 8 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i 2 i dτ ≤ 1 4 t 0 Γi,j ν i ∂ n n n i u k i -b b b i • n n n i u k i + S i,j u k i 2 dσ dτ + t 0 Γi,j (p i,j + p j,i )(-p i,j + p j,i + b b b i •n n n i )(u k i ) 2 dσ dτ + q 2 ( b b b i i,1,∞ + c i i,∞ ) u k i (t) 2 i + q 2 t 0 Γi,j (-p i,j +p j,i +2b b b i •n n n i )(∂ t u k i +∇ Γi,j •(r r r i,j u k i )-∇ Γi,j •(s i,j ∇ Γi,j u k i )) u k i dσ dτ + C t 0 β i u k i 2 i dτ + q t 0 √ ν i ∇u k i 2 i dτ . (2.28)
In order to estimate the fourth term in the right-hand side of (2.28), we observe that

t 0 Γi,j (-p i,j + p j,i + 2b b b i • n n n i )u k i ∂ t u k i dσ dτ = 1 2 Γi,j (-p i,j + p j,i + 2b b b i • n n n i )u k i (t) 2 dσ.
By the trace theorem in the right-hand side, we write:

t 0 Γi,j (-p i,j + p j,i + 2b b b i • n n n i )u k i ∂ t u k i dσ dτ ≤ C u k i (t) i √ ν i ∇u k i (t) i , and
u k i (t) 2 i = 2 t 0 Ωi (∂ t u k i )u k i ≤ 2 t 0 ∂ t u k i 2 i 1 2 t 0 u k i 2 i 1 2 , (2.29)
we obtain

q 2 t 0 Γi,j (-p i,j + p j,i + 2b b b i • n n n i )u k i ∂ t u k i dσ dτ ≤ q 8 t 0 ∂ t u k i 2 i dτ + q 4 u k i (t) 2 i +C t 0 β i u k i 2 i dτ . (2.30)
Moreover, integrating by parts and using the trace theorem, we have

- q 2 t 0 Γi,j ∇ Γi,j • (s i,j ∇ Γi,j u k i )(-p i,j + p j,i + 2b b b i • n n n i )u k i dσ dτ ≤ q 16 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i 2 i dτ + C( t 0 ∇ Γi,j u k i 2 i dτ + t 0 β i u k i 2 i dτ ). (2.31)
Using (2.29), we estimate the third term in the right-hand side of (2.28) by

q 2 ( b b b i i,1,∞ + c i i,∞ ) u k i (t) 2 i ≤ q 8 t 0 ∂ t u k i 2 i dτ + C t 0 β i u k i 2 i dτ. (2.32)
Replacing (2.31), (2.30) and (2.32) in (2.28), then using the transmission conditions, we have:

1 2 β i u k i (t) 2 i + q 2 u k i (t) 2 i +q ψ i,j √ s i,j ∇ Γi,j u k i (t) 2 i + t 0 β i u k i (τ, •) 2 i dτ + q 4 t 0 ∂ t u k i 2 i dτ + q 16 t 0 ψ i,j √ ν i s i,j ∇ ∇ Γi,j u k i 2 i dτ + 1 4 t 0 Γi,j ν i ∂ n n n i u k i -b b b i • n n n i u k i -S j,i u k i 2 dσ dτ ≤ 1 4 t 0 Γi,j ν j ∂ n n n i u k-1 j -b b b j • n n n i u k-1 j + S i,j u k-1 j 2 dσ dτ + C t 0 β i u k i 2 i dτ + q 2 t 0 √ ν i ∇u k i 2 i dτ .
We now sum up over the interfaces j ∈ N i , then over the subdomains for 1 ≤ i ≤ I, and on the iterations for 1 ≤ k ≤ K, the boundary terms cancel out, and we obtain for any t ∈ (0, T ),

k∈[[1,K]] i∈[[1,I]] β i,j u k i (t) 2 i + q √ ν i ∇u k i (t) 2 i + ν 0 t 0 ∇(β i,j u k i ) 2 i dτ ≤ α(t) + C k∈[[1,K]] i∈[[1,I]] t 0 β i,j u k i 2 i dτ + q t 0 √ ν i ∇u k i 2 i dτ , (2.33) with α(t) = 1 4 i∈[[1,I]] j∈Ni t 0 Γi,j ν j ∂ n n n i u 0 i -b b b j • n n n i u 0 j + S i,j u 0 j 2 dσ dτ. (2.34)
We conclude with Grönwall's lemma as before.

3. The discontinous Galerkin time stepping for the Schwarz waveform relaxation algorithm. In the following sections, in order to simplify the analysis, we suppose that c

+ 1 2 ∇ • b b b ≥ α 0 > 0 a.e. in Ω.
3.1. Time discretization of the local problem: discontinuous Galerkin method. We suppose that the coefficients are restricted to p -b b b•n n n 2 + q 2 ∇ Γ • r r r > 0 a.e. on Γ, q ≥ 0 a.e. and s > 0 a.e.. This implies that the bilinear form a defined in (2.6) is positive definite on H 1 (O) when q = 0, and positive definite on H 1 1 (O) when q ≥ q 0 > 0 a.e.

We recall the time-discontinuous Galerkin method, as presented in [START_REF] Johnson | Time discretization of parabolic problems by the discontinuous Galerkin method[END_REF]. We are given a decomposition T of the time interval (0, T ), I n = (t n , t n+1 ], for 0 ≤ n ≤ N , the mesh size is k n = t n+1t n . For B a Banach space and I an interval of R, define for any integer d ≥ 0

P d (B) = {ϕ : I → B, ϕ(t) = d i=0 ϕ i t i , ϕ i ∈ B}, P d (B, T ) = {ϕ : I → B, ϕ |In ∈ P d (B), 1 ≤ n ≤ N }. Let B = H 1 1 (O) if q > 0, B = H 1 (O) if q = 0.
We define an approximation U of u, polynomial of degree lower than d on every subinterval I n . For every point t n , we define U (t - n ) = lim t→tn-0 U (t), and note U (t + n ) = lim t→tn+0 U (t). The time discretization of (2.7) leads to searching

U ∈ P d (B, T ) such that            U (0, •) = u 0 , ∀ V ∈ P d (B, T ) : In ( m( U , V ) + a(U, V )) dt +m(U (t + n , •) -U (t - n , •), V (t + n , •)) = In L(V ) dt, (3.1)
with

L(V ) = (f, V ) L 2 (O) + (g, V ) L 2 (Γ) . Since I n is closed at t n+1 , U (t - n+1
) is the value of U at t n+1 . Due to the discontinuous nature of the test and trial spaces, the method is an implicit time stepping scheme, and U ∈ P d (B, T ) is obtained recursively on each subinterval, which makes it very flexible.

Theorem 3.1. If p -b b b•n n n 2 + q 2 ∇ Γ • r r r > 0 a.e. on Γ, q ≥ 0 a.e. and s > 0 a.e., equation (3.1) defines a unique solution.

Proof. The result relies on the fact that the bilinear form a is definite positive. It is is most easily seen by using a basis of Legendre polynomials. U ∈ P d (H 1 1 (Ω), T ) is obtained recursively on each subinterval. We introduce the Legendre polynomials L n , orthogonal basis in L 2 (-1, 1), with L n (1) = 1. L n has the parity of n, hence L n (-1) = (-1) n . A basis of orthogonal polynomial on I n is given by

L n,k (t) = L k ( 2 kn (t -tn+1+tn 2
)). Choose V (t, x) = L n,j (t)Φ j (x) in (3.1) with Φ j ∈ H 1 1 (Ω), and expand U on I n as U (t, x) = d k=0 U k (x)L n,k (t). Suppose U to be given on (0, t n ].
In order to determine U on I n , we must solve the system: for any Φ j ∈ H 1 1 (Ω),

d k=0 In Ln,k L n,j m(U k , Φ j ) + L n,k L n,j a(U k , Φ j ) dt + d k=0 L n,k (t + n )L n,j (t + n )m(U k , Φ j ) = In L n,j L(Φ j ) dt.
It is an implicit scheme. We calculate the coefficients

In L n,k L n,j = δ kj L n,j 2 , In Ln,k L n,j = 0 if k ≤ j 1 -(-1) k+j if k > j, In Ln,k L n,j + L n,k (t + n )L n,j (t + n ) = (-1) k+j if k ≤ j 1 if k > j,
which leads to

Ln,j 2 a(Uj, Φj ) + m(Uj , Φj ) + k<j (-1) k+j m(U k , Φj ) + k>j m(U k , Φj ) = In Ln,j L(Φj ) dt.
It is a square system of partial differential equations, of the type coercive + compact. By the Fredholm alternative, we only need to prove uniqueness. Choose now Φ j = U j , and obtain

j L n,j 2 a(U j , U j ) + j m(U j , U j ) + 2 j k>j k+j even m(U k , U j ) = 0,
and since a is positive definite, we deduce that U = 0. We will make use of the following remark ( [START_REF] Makridakis | A posteriori error analysis for higher order dissipative methods for evolution problems[END_REF]). We introduce the Gauss-Radau points, (0 < τ 1 , . . . , τ d+1 = 1), defined such that the quadrature formula

1 0 f (t)dt ≈ d+1 j=1 w q f (τ q )
is exact in P 2d , and the interpolation operator I n on [t n , t n+1 ] at points (t n , t n + τ 1 k n , . . . , t n + τ d+1 k n ). For any χ ∈ P d , χ = I n χ ∈ P d+1 .

Let I : P d (B, T ) → P d+1 (B, T ) be the operator whose restriction to each subinterval is I n and satisfies IU (t + n ) = U (t - n ). By using the Gauss-Radau formula, which is exact in P 2d , we have for all ψ i,j ∈ P d In> dIχ dt ψ i,j dt -

In dχ dt ψ i,j dt = (χ(t + n ) -χ(t - n ))ψ i,j (t + n ).
As a consequence, we have a very useful inequality:

In d dt (Iψ i,j )ψ i,j dt ≥ 1 2 [ψ i,j (t - n+1 ) 2 -ψ i,j (t - n ) 2 ]. (3.2)
Also, equation (3.1) can be rewritten as

In ( m( dIU dt , V ) + a(U, V )) dt = In L(V ) dt, (3.3)
or in the strong formulation:

∂ t (IU ) + ∇ • (b b bU -ν∇U ) + cU = P f, in Ω × (0, T ), ν ∂ n n n -b b b • n n n U + p U + q(∂ t (IU ) + ∇ Γ • (r r rU -s∇ Γ U )) = P g on Γ × (0, T ). (3.4)
Here P is the projection L 2 in each subinterval of T on P d .

Theorem 3.2 (Thomee, [START_REF] Thomee | Galerkin Finite Element Methods for Parabolic Problems[END_REF]). Let U be the solution of (3.1) and u the solution of (2.1). Under the assumptions of Theorem 2.1, the estimate holds

u -U L ∞ (In,L 2 (Ω)) ≤ Ck d+1 ∂ d+1 t u L 2 (0,T ;H 2 2 (Ω)) , (3.5
)

with k = max 0≤n≤N k n .
3.2. The discrete in time optimized Schwarz waveform relaxation algorithm with different subdomains grids. In this part we present and analyse the discrete non conforming in time optimized Schwarz waveform relaxation method.

The time partition in subdomain Ω i , is T i , with N i + 1 intervals I i n , and mesh size k i n . In view of formulation (3.4), we define interpolation operators I i and projection operators P i in each subdomain, i.e. P i is the projection L 2 in each subinterval of T i on P d , and we solve

∂ t (I i U k i ) + ∇ • (b b b i U k i -ν i ∇U k i ) + c i U k i = P i f in Ω i × (0, T ), (3.6a
)

ν i ∂ n n n i -b b b i • n n n i U k i + S i,j U k i = P i (ν j ∂ n n n i -b b b j • n n n i ) U k-1 j + S i,j U k-1 j on Γ i,j , j ∈ N i . (3.6b)
Here the operators are different on either part of the "equal" sign:

S i,j U = p i,j U + q i,j (∂ t (I i U ) + ∇ Γi,j • (r r r i,j U -s i,j ∇ Γi,j U )) S i,j U = p i,j U + q i,j (∂ t (I j U ) + ∇ Γi,j • (r r r i,j U -s i,j ∇ Γi,j U )). (3.7)
Formally, the sequence of problems (3.6) converges to the solution of

∂ t (I i U i ) + ∇ • (b b b i U i -ν i ∇U i ) + c i U i = P i f in Ω i × (0, T ), (3.8a
)

ν i ∂ n n n i -b b b i • n n n i U i + S i,j U i = P i (ν j ∂ n n n i -b b b j • n n n i ) U j + S i,j U j on Γ i,j , j ∈ N i . (3.8b)
We present the analysis first with Robin transmission conditions (e.g. q i,j = 0) and general decomposition, and then with order 2 transmission conditions and decomposition in strips.

The Robin case.

We consider here a general decomposition of the domain, possibly with corners. We solve (3.6) with q i,j = 0, i.e. S i,j U = S i,j U = p i,j U .

Theorem 3.3. Assume q i,j = 0, p j,i -p i,j -b b b i • n n n i = 0, p i,j -b b bi•n n ni 2 > 0. Problem (3.8
) has a unique solution (U i ) i∈I , and U i is the limit of the iterates of algorithm (3.6).

Proof. We first write energy estimates on (3.6) for f ≡ 0 and u 0 ≡ 0. We start like in the proof of Theorem 2.2. We multiply (3.6a) by U k i , integrate on Ω i , then integrate on the interval (t i n , t i n+1 ) and use (3.2) and (2.14):

U k i (t i n+1 ) 2 L 2 (Ωi) -U k i (t i n ) 2 L 2 (Ωi)
+ 2

I i n (ν i ∇U k i , ∇U k i ) L 2 (Ωi) + ((c i + 1 2 ∇ • b b b i )U k i , U k i ) L 2 (Ωi) dτ + j∈Ni I i n Γi,j 1 p i,j + p j,i ν i ∂ n n ni u k i -b b b i • n n n i U k i -p j,i U k i 2 dσ dτ ≤ j∈Ni I i n Γi,j 1 p i,j + p j,i ν i ∂ n n ni u k i -b b b i • n n n i U k i + p i,j U k i 2 dσ dτ + j∈Ni I i n Γi,j (p j,i -p i,j -b b b i • n n n i )(U k i ) 2 dσ dτ.
We can not use Grönwall's Lemma like in the continuous case, due to the presence of the global in time projection operator P j in the transmission conditions. Therefore we have to assume that p j,ip i,j -b b b i •n n n i = 0 everywhere, which cancels the last term. We sum up over the time intervals, using the fact that the errors vanish at time 0:

U k i (T ) 2 L 2 (Ωi) + 2 min(ν 0 , α 0 ) T 0 U k i 2 H 1 (Ωi) dτ + j∈Ni T 0 Γi,j 1 p i,j + p j,i ν i ∂ n n n i u k i -b b b i • n n n i U k i -p j,i U k i 2 dσ dτ ≤ j∈Ni T 0 Γi,j 1 p i,j + p j,i ν i ∂ n n n i u k i -b b b i • n n n i U k i + p i,j U k i 2 dσ dτ.
We now insert the transmission conditions

U k i (T ) 2 L 2 (Ωi) + 2 min(ν 0 , α 0 ) T 0 U k i 2 H 1 (Ωi) dτ + j∈Ni T 0 Γi,j 1 p i,j + p j,i ν i ∂ n n n i u k i -b b b i • n n n i U k i -p j,i U k i 2 dσ dτ ≤ j∈Ni T 0 Γi,j 1 p i,j + p j,i P i ν j ∂ n n ni U k-1 j -b b b j • n n n i U k-1 j + p i,j U k-1 j 2 dσ dτ.
We use the fact that the projection operator is a contraction to obtain:

U k i (T ) 2 L 2 (Ωi) + 2 min(ν 0 , α 0 ) T 0 U k i 2 H 1 (Ωi) dτ + j∈Ni T 0 Γi,j 1 p i,j + p j,i ν i ∂ n n n i u k i -b b b i • n n n i U k i -p j,i U k i 2 dσ dτ ≤ j∈Ni T 0 Γj,i 1 p i,j + p j,i ν j ∂ n n n j U k-1 j -b b b j • n n n j U k-1 j -p i,j U k-1 j 2 dσ dτ.
We sum up over the subdomains, we define the boundary term

B k = i∈[[1,I]] j∈Ni T 0 Γi,j 1 p i,j + p j,i ν i ∂ n n ni u k i -b b b i • n n n i U k i -p j,i U k i 2 dσ dτ, we obtain i∈[[1,I]] ( U k i (T ) 2 L 2 (Ωi) + 2 min(ν 0 , α 0 ) T 0 U k i 2 H 1 (Ωi) dτ ) + B k ≤ B k-1 . (3.9)
We first apply this inequality to prove the first part of the Theorem. (3.8) is a square discrete system, and proving well-posedness is equivalent to proving uniqueness. Dropping the superscript in (3.9) gives the result. As for the convergence, we proceed as in the continuous case by summing (3.9) over the iterates to obtain

that i∈[[1,I]] U k i (T ) 2 L 2 (Ωi) and i∈[[1,I]] T 0 U k i 2
H 1 (Ωi) dτ tend to zero as k tend to infinity.

3.2.2. The Order 2 case. We restrict ourselves to a splitting of the domain into strips with parallel planar interfaces.

Theorem 3.4. We assume that p i,j = p > 0, q i,j = q > 0, s i,j = s > 0, b b b i = 0 and r r r i,j = 0. Problem (3.8) has a unique solution (U i ) i∈I , and U i is the limit of the iterates of algorithm (3.6).

Proof. We consider the algorithm (3.6) on the error, so we suppose f = u 0 = 0. As in the continuous case, the proof is based on energy estimates containing the term

I i n Γi,j ν i ∂ n n n i u k i + S i,j U k i 2 dσ dτ,
and that we derive by multiplying successively the first equation of (3.6) by the terms

U k i , ∂ t (I i U k i ), and -∆ Γi,j U k i . We set ϕ 2 i = √ ν i ∇ϕ 2 L 2 (Ωi) + √ c i ϕ 2 L 2 (Ωi)
. We multiply the first equation of (3.6) by U k i , we integrate on I i n × Ω i then integrate by parts in space and use (3.2):

1 2 U k i (t - n+1 ) 2 i + I i n U k i (τ, •) 2 i dτ - I i n Γi,j ν i ∂ n n ni u k i U k i dσ dτ ≤ 1 2 U k i (t - n ) 2 i .
(3.10)

We multiply the first equation of (3.6) by ∂ t (I i U k i ), integrate on I i n × Ω i and then integrate by parts in space and use (3.2):

1 2 U k i (t - n+1 ) 2 i + I i n ∂ t (I i U k i ) 2 i dτ - I i n Γi,j ν i ∂ n n ni u k i ∂ t (I i U k i ) dσ dτ ≤ 1 2 U k i (t - n ) 2 i . (3.11)
Now we multiply the first equation of (3.6) by -∆ Γi,j U k i integrate on I i n × Ω i and integrate by parts in space and use (3.2):

1 2 ∇ Γi,j U k i (t - n+1 ) 2 i + I i n ∇ Γi,j U k i (t, .) 2 i + I i n Γi,j ν i ∂ n n n i u k i ∆ Γi,j u k i dσ dτ ≤ 1 2 ∇ Γi,j U k i (t - n ) 2 i , (3.12)
where we have used the fact that ∆ Γi,j is a constant coefficient operator. Let

E n (U k i ) = p 2 U k i ((t i n ) -) 2 i + q 2 U k i ((t i n ) -) 2 i + sq 2 ∇ Γi,j U k i ((t i n ) -) 2 i .
Multiplying (3.10) by p, (3.11) by q and (3.12) by sq, and adding the three equations with (3.10), we get

E n+1 (U k i ) + I i n [ p U k i (t, •) 2 i +q ∂ t (I i n U k i ) 2 i + sq ∇ Γi,j U k i (t, •) 2 i ] dt - j∈Ni I i n Γij ν i ∂ n n n i u k i S ij U k i dx 2 dt ≤ E n (U k i ).
It can be rewritten as

E n+1 (U k i ) + I i n [ U k i (t, •) 2 i +q ∂ t (I i n U k i ) 2 i + sq ∇ Γi,j U k i (t, •) 2 i ] dt + 1 4 j∈Ni I i n Γij (ν i ∂ n n n i u k i -S ij U k i ) 2 ≤ E n (U k i )+ 1 4 j∈Ni I i n Γij (ν i ∂ n n n i u k i + S ij U k i ) 2 .
We now sum in time for 0 ≤ n ≤ N , and use the transmission condition. Since E 0 (U k i ) = 0, we obtain

E N +1 (U k i ) + T 0 [ U k i (t, •) 2 i +q ∂ t (I i n U k i ) 2 i + sq ∇ Γi,j U k i (t, •) 2 i ] dt + 1 4 j∈Ni T 0 Γij (ν i ∂ n n n i u k i -S ij U k i ) 2 dt ≤ 1 4 j∈Ni T 0 Γij (P i (-ν j ∂ n n n j U k-1 j + S ij U k-1 j )) 2 dt.
We sum up over the subdomains and use the fact that the projection is a contraction. Since we are in the case where p ij = p, q ij = q, r ij = 0 and s ij = s, we have S ij = S ji . Thus, we can sum up over the iterates, the boundary terms cancel out, and we obtain

K k=1 I i=1 E N +1 (U k i ) + T 0 [ U k i (t, •) 2 i +q ∂ t (I i n U k i ) 2 i + sq ∇ Γi,j U k i (t, •) 2 i ] dt + 1 4 I i=1 T 0 Γi (ν i ∂ n n n i u k i -S ij U K i ) 2 ≤ 1 4 I i=1 T 0 Γi (ν i ∂ n n n i U 0 i -S ij U 0 i ) 2 .
We conclude as in the proof of Theorem 3.4. We now state the error estimate in the Robin case.

3.3.

Error estimates in the Robin case. Theorem 3.5.

If ∇ • b b b = 0, p i,j -b b bi•n n ni 2 = p j,i - b b bj •n n nj 2
= p > 0, and q i,j = 0, the error between u and the solution U i of (3.8) is estimated by:

I i=1 u -U i 2 L ∞ (0,T,L 2 (Ωi)) ≤ Ck 2(d+1) ∂ d+1 t u 2 L 2 (0,T ;H 2 (Ω)) . (3.13)
Proof. We introduce the projection operator

P - i as ∀n ∈ [[1, N i ]], P - i ϕ ∈ P d (I i n ), P - i ϕ(t i n+1 ) = ϕ(t i n+1 ), ∀ψ i,j ∈ P d-1 (I i n ), I i n (P - i ϕ -ϕ)(t)ψ i,j (t) dt = 0.
We define W i = P - i (u| Ωi ), Θ i = U i -W i and ρ i = W i -u| Ωi . Classical projection estimates ( [START_REF] Thomee | Galerkin Finite Element Methods for Parabolic Problems[END_REF]) yield the estimate on ρ i :

I i=1 ρ i 2 L ∞ (0,T,L 2 (Ωi)) ≤ Ck 2(d+1) ∂ d+1 t u 2 L 2 (0,T ;L 2 (Ω)) .
Thus, since U i -u| Ωi = Θ i + ρ i , it suffices to prove estimate (3.13) for Θ i . Now, thanks to the equations on u and U i , and the identity d dt

I i P - i = P i d dt , Θ i satisfies: ∂ t (I i Θ i ) + b b b • ∇Θ i -ν∆Θ i + cΘ i = -b b b • ∇ρ i + ν∆ρ i -cρ i +(1 -P i )(∂ t u -f ) in Ω i × (0, T ), ν i ∂ n n n i -b b b i • n n n i Θ i + p i,j Θ i = P i ( ν j ∂ n n n i -b b b j • n n n i Θ j + p i,j Θ j ) -(1 -P i )( ν j ∂ n n ni -b b b j • n n n i W j + p i,j W j ) on Γ ij × (0, T ), j ∈ N i . (3.14)
Multiply the first equation of (3.14) by Θ i , integrate on (t i n , t i n+1 ) × Ω i , using (3.2) and integrate by parts in space. Terminate with Cauchy Schwarz inequality:

1 2 Θ i ((t i n+1 ) -) 2 i + I i n Θ i (t, •) 2 i dt - I i n Γi (ν i ∂ n n n i Θ i - b b b i • n n n i 2 Θ i ) Θ i dx 2 dt ≤ 1 2 Θ i ((t i n ) -) 2 i + C I i n ρ i (t, •) 2 H 2 (Ωi) dt.
Rewriting the boundary integral using (2.14), we obtain

1 2 Θ i ((t i n+1 ) -) 2 i + I i n Θ i (t, •) 2 i dt+ 1 4p j∈Ni I i n Γij (ν i ∂ n n ni Θ i -b b b i •n n n i Θ i -p j,i Θ i ) 2 dx 2 dt ≤ 1 4p j∈Ni I i n Γij (ν i ∂ n n ni Θ i -b b b i •n n n i Θ i +p i,j Θ i ) 2 dx 2 + 1 2 Θ i ((t i n ) -) 2 i +C I i n ρ i (t, •) 2 H 2 (Ωi) dt.
Using the transmission condition in (3.14) together with the fact that P i and 1 -P i are orthogonal to each other and have norm 1, we get by a trace theorem

1 2 Θ i ((t i n+1 ) -) 2 i + I i n Θ i (t, •) 2 i dt+ 1 4p j∈Ni I i n Γij (ν i ∂ n n n i Θ i -b b b i •n n n i Θ i -p j,i Θ i ) 2 dx 2 dt ≤ 1 4p j∈Ni I i n Γij (ν j ∂ n n n j Θ j -b b b j • n n n j Θ j -p i,j Θ i ) 2 dx 2 + 1 2 Θ i ((t i n ) -) 2 i + C I i n ρ i (t, •) 2 H 2 (Ωi) dt + C I i n (1 -P i )(u| Ωi )(t, •) 2 H 2 (Ωi) dt. (3.15)
Classical error estimates [START_REF] Thomee | Galerkin Finite Element Methods for Parabolic Problems[END_REF] imply:

T 0 ρ i (t, •) 2 H 2 (Ωi) dt+ T 0 (1-P i )(u| Ωi )(t, •) 2 H 2 (Ωi) ≤ Ck 2(d+1) ∂ d+1 t u 2 L 2 (0,T ;H 2 (Ωi)) .
Summing (3.15) in i and n, and using the previous equation yields (3.13).

4. Space-time nonconforming algorithm. In this section we describe the implementation of algorithm (3.6), especially in the cases d = 0 and d = 1. We start from the semi-disrete in time scheme and use finite elements for the space discretization in each subdomain. In order to permit non-matching grids in time and space on the boundary, we extend the nonconforming approach in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF].

We describe first the implementation of algorithm (3.6) at the semi-discrete in time level, and then at the space-time discret level.

4.1. Time discretization. We recall the subdomain scheme in time, and give it in details for d = 0 and d = 1. Then we describe the computation of the transmission conditions in algorithm (3.6).

4.1.1. Interior scheme. We consider the subdomain problem in the algorithm (3.6)

at iteration k in O = Ω i . Let B i = H 1 1 (Ω i ) if q > 0, B i = H 1 (Ω i ) if q = 0. We set U = U k i ∈ P d (B i , T i
), and we omit the subscript i for the local time scheme to simplify the notations :

∂ t (IU ) + ∇ • (b b bU -ν∇U ) + cU = P f in O × (0, T ) ν ∂ n n n -b b b • n n n U + p U + q(∂ t (IU ) + ∇ Γ • (r r rU -s∇ Γ U )) = P g on Γ × (0, T ). ( 4

.1)

Case d = 0

In the case d = 0, the approximating functions are piecewise constant in time, then

U (t) = U n+1 = U n + in I n , we have I n U = U n + t-tn kn (U n+1 -U n ), P ξ = 1
kn In ξ(•, s) ds and the method reduces to the modifed backward Euler method

U n+1 k n + ∇ • (b b bU n+1 -ν∇U n+1 ) + cU n+1 = U n k n + 1 k n In f (•, s) ds in O ν ∂ n n n -b b b • n n n U n+1 + p U n+1 + q( U n+1 k n + ∇ Γ • (r r rU n+1 -s∇ Γ U n+1 )) = q U n k n + 1 k n In g(•, s) ds on Γ. (4.2) Case d = 1
In that case, for piecewise linear functions of t, using a basis of Legendre polynomials we may write,

U (t) = U n+1 0 + 2 t-t n+1/2 kn U n+1 1 , on I i n , with t n+1/2 = tn+tn+1 2 , U n = U (t n
), and we have on I i n :

I n U = 1 4 (5U n+1 0 -U n+1 1 -U n ) + (U n+1 0 + U n+1 1 -U n ) t -t n+1/2 k n +3(-U n+1 0 + U n+1 1 + U n ) t -t n+1/2 k n 2 ,
and

P ξ = ξ 0 + 2 t-t n+1/2 kn ξ 1 with ξ 0 = 1 kn In ξ(•, s) ds ξ 1 = 6 kn In s-t n+1/2 kn ξ(•, s) ds (4.3)
Thus, we obtain for the determination of U n+1 0 and U n+1

1 the system 1 k n (U n+1 0 + U n+1 1 ) + ∇ • (b b bU n+1 0 -ν∇U n+1 0 ) + cU n+1 0 = U n k n + 1 k n In f (•, s) ds 3 k n (-U n+1 0 + U n+1 1 ) + ∇ • (b b bU n+1 1 -ν∇U n+1 1 ) + cU n+1 1 = - 3U n k n + 6 k n In s -t n+1/2 k n f (•, s) ds in O ν ∂ n n n -b b b • n n n U n+1 0 + p U n+1 0 + q( 1 k n (U n+1 0 + U n+1 1 ) + ∇ Γ • (r r rU n+1 0 -s∇ Γ U n+1 0 )) = q U n k n + 1 k n In g(•, s) ds ν ∂ n n n -b b b • n n n U n+1 1 + p U n+1 1 + q( 3 k n (-U n+1 0 + U n+1 1 ) + ∇ Γ • (r r rU n+1 1 -s∇ Γ U n+1 1 )) = -q 3U n k n + 6 k n In s -t n+1/2 k n g(•, s) ds on Γ.
(4.4) Multiplying the first equation of (4.2) by v ∈ B i (resp. the first equation of (4.4) by v ∈ B i and the second equation of (4.4) by w ∈ B i ), integrating by parts on O, and using the boundary conditions, the variational formulation is:

Case d = 0 (Variational formulation) m(U n+1 , v) + k n a(U n+1 , v) = m(U n , v) + I i n (f (•, s), v)ds + I i n (g(•, s), v) Γ ds, ∀v ∈ B i . (4.5) Case d = 1 (Variational formulation) m(U n+1 0 , v) + k n a(U n+1 0 , v) + m(U n+1 1 , v) = m(U n , v) + In (f (•, s), v)ds + In (g(•, s), v) Γ ds, -m(U n+1 0 , v) + m(U n+1 1 , w) + k n a(U n+1 1 , w) = -m(U n , v) + In 2(s -t n+1/2 ) k n (f (•, s), w)ds + In 2(s -t n+1/2 ) k n (g(•, s), w) Γ ds, ∀v, w ∈ B i . (4.6)
Remark 4.1. Equations (4.5) and (4.6) can be derived directly from (3.1). However we will need formulas (4.2) and (4.4) in the space nonconforming case.

We now discuss the computation of the right-hand side on the interface Γ i,j ×(0, T ) for j ∈ N i in the algorithm (3.6).

Transmission terms. Let (g 1

i,j ) be a given initial guess in P d (L 2 (Γ i,j ), T i ), for 1 ≤ i ≤ I, j ∈ N i . Then, at iteration k ≥ 1, we solve the subdomain problem in Ω i :

∂ t (I i U k i ) + ∇ • (b b b i U k i -ν i ∇U k i ) + c i U k i = P i f in Ω i × (0, T ), (4.7a)
ν i ∂ n n n i -b b b i • n n n i U k i + S i,j U k i = g k i,j , j ∈ N i on Γ i,j × (0, T ), (4.7b)
The function g k i,j is defined for k ≥ 2 by

g k i,j = P i gk j,i , (4.8)
with gk j,i , k ≥ 2, defined by

gk j,i = -(ν j ∂ n n n j -b b b j • n n n j ) U k-1 j + S i,j U k-1 j .
We remark that, for k ≥ 2,

gk j,i = -g k-1 j,i + S j,i U k-1 j + S i,j U k-1 j = -g k-1 j,i + (p i,j + p j,i )U k-1 j + (q i,j + q j,i )∂ t (I j U k-1 j ) + q i,j (∇ Γi,j • (r r r i,j U k-1 j -s i,j ∇ Γi,j U k-1 j )) + q j,i (∇ Γj,i • (r r r j,i U k-1 j -s j,i ∇ Γj,i U k-1 j )).
Once gk j,i is computed from U k-1 j , we obtain g k i,j from (4.8) as follows : we introduce the basis functions (ϕ i n,α ) 0≤α≤d of polynomial of degree lower than d on subinterval

I i n , then (g k i,j ) |I i n = (P i gk j,i ) |I i n = d α=0 G i,k n,α ϕ i n,α with G i,k n,α ∈ L 2 (Γ i,j ) solution of the system d α=0 G i,k n,α I i n ϕ i n,α ϕ i n,β ds = I i n gk j,i ϕ i n,β , β ∈ {0, ..., d}.
Thus, the computation of g k i,j on each I i n needs the computation of terms in the form

I i n gk j,i ϕ i n,β ds, (4.9)
for β ∈ {0, ..., d}. Recall that gk j,i is defined on Γ i,j × I and gk j,i ∈ P d (L 2 (Γ i,j ), T j ). Thus, we first write the integral in (4.9) as an integral over I : let Φ i n,α be the function defined on I, equal to ϕ i n,α on I i n and equal to zero on I\I i n . Then

I i n gk j,i ϕ i n,β ds = I gk j,i Φ i n,β ds. (4.10)
We now decompose gk j,i on the basis functions (ϕ j m,α ) 0≤α≤d of polynomial of degree lower than d on each subinterval I j m :

(g k j,i ) |I j m = d α=0 Gj,k m,α ϕ j m,α , with Gj,k m,α ∈ L 2 (Γ i,j ) solution of the system d α=0 Gj,k m,α I j m ϕ j m,α ϕ j m,β ds = I j m gk j,i ϕ j m,β , β ∈ {0, ..., d}.
Introducing the function Φ j m,α defined on I, equal to ϕ j m,α on I j m and equal to zero on I\I j m , we have

gk j,i = Nj m=0 d α=0 Gj,k m,α Φ j m,α , (4.11)
Replacing (4.11) in (4.10) leads to

I i n gk j,i ϕ i n,β ds = Nj m=0 d α=0 Gj,k m,α I Φ j m,α Φ i n,β ds.
Let M α,β be the projection matrix defined by

(M α,β ) n+1,m+1 = I Φ j m,α Φ i n,β ds. 0 ≤, n ≤ N i , 0 ≤ m ≤ N j .
Then we have, for 0 ≤ n ≤ N i ,

I i n gk j,i ϕ i n,β ds = d α=0 (M α,β Gj,k α ) n with Gj,k α = ( Gj,k 0,α , ..., Gj,k Nj ,α ) t .
In the special cases d = 0 and d = 1, we obtain :

d = 0
In that case there is one basis function ϕ i n,0 = 1 on I i n , and

I i n gk j,i ϕ i n,0 ds = I i n gk j,i ds = (M 0,0 Gj,k 0 ) n , with (M 0,0 ) n+1,m+1 = I 1 I j m 1 I i n ds, Gj,k m,0 = 1 k j m I j m gk j,i ds, 0 ≤, n ≤ N i , 0 ≤ m ≤ N j . Case d = 1
In that case there are two basis functions

ϕ i n,0 = 1, ϕ i n,1 = 2 s-t i n+1/2 k i n on I i n ,
and

I i n gk j,i ϕ i n,0 ds = (M 0,0 Gj,k 0 + M 1,0 Gj,k 1 ) n , I i n gk j,i ϕ i n,1 ds = (M 0,1 Gj,k 0 + M 1,1 Gj,k 1 ) n , with, for 0 ≤, n ≤ N i , 0 ≤ m ≤ N j , (M 0,0 ) n+1,m+1 = I 1 I j m 1 I i n ds, (M 1,1 ) n+1,m+1 = 4 I s -t j m+1/2 k j m 1 I j m s -t i n+1/2 k i n 1 I i n ds, (M 1,0 ) n+1,m+1 = 2 I s -t j m+1/2 k j m 1 I j m 1 I i n ds, (M 0,1 ) n+1,m+1 = 2 I 1 I j m s -t i n+1/2 k i n 1 I i n ds,
and Gj,k m,0 , Gj,k m,1 defined by

   Gj,k m,0 = 1 k j m I j m gk j,i ds Gj,k m,1 = 6 k j m I j m s-t j m+1/2 k j m gk j,i ds
The projection matrices M α,β are computed by a simple and optimal projection algorithm without any additional grid (see [START_REF] Gander | Optimal Schwarz waveform relaxation for the one dimensional wave equation[END_REF], [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF]).

We now discuss the space dicretization using finite elements.

4.2. Space discretization. We suppose that each subdomain Ω i is provided with its own mesh T i h , 1 ≤ i ≤ I, such that Ω i = ∪ T ∈T i h T. For T ∈ T i h , let h T := sup x,y∈T d(x, y) be the diameter of T and h the discretization parameter

h = max 1≤i≤I h i , with h i = max T ∈T i h h T .
Let P 1 (T ) denote the space of all polynomials defined over T of total degree less than or equal to 1. Then, we define over each subdomain the conforming spaces V i h by :

V i h = {v i,h ∈ C 0 (Ω i ), v i,h |T ∈ P 1 (T ), ∀T ∈ T i h }.
In what follows we assume that the mesh is designed by taking into account the geometry of the Γ i,j in the sense that, the space of traces over each Γ i,j of elements of V i h is a finite element space denoted by V i,j h . Let n i,j be the dimension of V i,j h and (χ i,j ℓ,h ) 1≤ℓ≤n i,j the finite element basis functions of V i,j h . We consider two cases : when the grids in space are conforming, and the case of nonconforming space grids. 4.2.1. Conforming case. In the case of conforming grids in space, we have V i,j h = V j,i h . We can replace B i by V i h in the variational formulation. We set :

ãi (u, v) = Ωi ( 1 2 ((b b b i • ∇u)v -(b b b i • ∇v)u)) dx + Ωi ν i ∇u • ∇v dx + Ωi (c i + 1 2 ∇ •b b b i )uv dx, (4.12
) and

< C i,j u, v > Γi,j = Γi,j (p i,j -b b bi•n n ni 2 ) uv +q i,j (∂ t (I i u) + ∇ Γi,j • (r r r i,j u))v -s i,j ∇ Γi,j u∇ Γi,j v dσ, < C i,j u, v > Γi,j = Γi,j (p i,j -b b bi•n n ni
2 ) uv +q i,j (∂ t (I j u) + ∇ Γi,j • (r r r i,j u))vs i,j ∇ Γi,j u∇ Γi,j v dσ.

We introduce the discret algorithm : let (g 1 i,j,h ) be a given initial guess in P d (V i,j h , T i ), for 1 ≤ i ≤ I, j ∈ N i . Let U k i,h be the approximation of u k i in P d (V i h , T j ). Then, at iteration k ≥ 1, we solve the subdomain problem in Ω i :

Ωi ∂ t (I i U k i,h )v i,h + ãi (U k i,h , v i,h ) dx + < C i,j U k i,h , v i,h > Γi,j = Ωi P i f v i,h dx + Γi,j g k i,j,h v i,h dσ, in (0, T ), ∀v i,h ∈ V i h , (4.13) For k ≥ 2, v h ∈ V i,j h , we define Γi,j g k i,j,h v h dσ := P i Γi,j gk j,i,h v h dσ, (4.14)
with Γi,j

gk j,i,h v h dσ := - Γi,j g k-1 j,i,h v h dσ+ < C j,i U k-1 j,h + C i,j U k-1 j,h , v h > Γi,j .
In equation (4.14) we used the fact that the space of traces over each Γ i,j of elements of V i h is the same as the space of traces over each Γ i,j of elements of V j h . For the computation of the right-hand side in (4.13), we follow the same steps as in section 4.1.2, where we replace gk j,i with g g g k j,i,h ∈ V The discrete formulation in the cases d = 0 and d = 1 are obtained from (4.5) and (4.6), by replacing B i by V i h .

When the space grids are nonconforming, following [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF], we cannot replace directly B i by the finite element space V i h in the variational formulation. We have to consider equation (4.1) (i.e. (4.2) for d = 0, and (4.4) for d = 1).

Nonconforming case.

In this section we extend the nonconforming approach in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF]. We consider the mortar spaces W i,j h as in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF]. Let m i,j be the dimension of W i,j h and (ψ i,j k,h ) 1≤k≤m i,j the finite element basis functions of W i,j h . We introduce the discrete algorithm : let

(U k-1 i,h , Q k-1 i,h ) ∈ P d (V i h , T i ) × P d ( W i,j h , T i ) be a discrete approximation of (U k-1 i , ν i ∂ n n n i U k-1 i ) in Ω i at step k -1. Then (U k i,h , Q k i,h ) is the solution in P d (V i h , T i ) × P d ( W i,j h , T i ) of d dt (I i U k i,h , v i,h ) i + ãi (U k i,h , v i,h ) i + Γi,j (Q k i,h - b b b i • n n n i 2 U k i,h )v i,h dσ = (P i f, v i,h ) i , in (0, T ), ∀v i,h ∈ V i h , Γi,j Q k i,h -b b b i • n n n i U k i,h + p i,j U k i,h ψ i,j h dσ + Γi,j q i,j (∂ t (I i U k i,h ) + ∇ Γi,j • (r r r i,j U k i,h ))ψ i,j h + q i,j s i,j ∇ Γi,j U k i,h ∇ Γi,j ψ i,j h dσ = Γi,j P i -Q k-1 j,h -b b b j • n n n i U k-1 j,h + p i,j U k-1 j,h ψ i,j h dσ + Γi,j P i q i,j (∂ t (I j U k-1 j,h ) + ∇ Γi,j • (r r r i,j U k-1 j,h ))ψ i,j h + q i,j s i,j ∇ Γi,j U k-1 j,h ∇ Γi,j ψ i,j h dσ on (0, T ), ∀ψ i,j h ∈ W i,j h , j ∈ N i . (4.16)
We give first the interior scheme for d = 0 and d = 1 and then the computation of the right-hand side in the transmission condition of (4.16).

Interior scheme. The discrete problem in subdomain O = Ω i in (4.16) is defined as follows : find

(U h , Q h ) := (U k i,h , Q k i,h ) in P d (V i h , T i ) × P d ( W i,j h , T i ) solution of d dt (IU h , v h ) + ã(U h , v h ) + Γ (Q h - b b b • n n n 2 U h )v h dσ = (P f, v h ), in (0, T ), ∀v h ∈ V i h , Γ (Q -b b b • n n nU h + p U h + q(∂ t (IU h ) + ∇ Γ • (r r rU h )))ψ h + qs∇ Γ U h ∇ Γ ψ h dσ = Γ (P g)ψ h dσ, on (0, T ), ∀ψ h ∈ W i,j h .
In the cases d = 0 and d = 1 we obtain

Case d = 0
In that case, the approximating functions are piecewise constant in time:

U h (t) = U n+1 h = U n h,+ and Q h (t) = Q n+1 h = Q n h,+ on I i n , and the discrete problem reduces to find (U n+1 h , Q n+1 h ) ∈ V i h × W i,j h solution of ( U n+1 h k n , v h ) + ã(U n+1 h , v h ) + Γ (Q n+1 h - b b b • n n n 2 U n+1 h )v h dσ = ( U n h k n , v h ) + 1 k n In (f (•, s), v h ) ds, ∀v h ∈ V i h , Γ Q n+1 h -b b b • n n nU n+1 h + p U n+1 h + q( U n+1 h k n + ∇ Γ • (r r rU n+1 h ))ψ h + qs∇ Γ U n+1 h ∇ Γ ψ h dσ = Γ q U n h k n ψ h dσ + 1 k n In Γ g(s)ψ h dσ ds, ∀ψ h ∈ W i,j h . Case d = 1
In that case, we write

U h (t) = U n+1 0,h + 2 t-t n+1/2 kn U n+1 1,h and Q h (t) = Q n+1 0,h + 2 t-t n+1/2 kn Q n+1 1,h on I i n , U n h = U h (t n ),
and the discrete problem reduces to find (U n+1 0,h , Q n+1 0,h ) and (U n+1

1,h , Q n+1 1,h ) in V i h × W i,j h solution of the system 1 k n (U n+1 0,h + U n+1 1,h , v h ) + ã(U n+1 0,h , v h ) + Γ (Q n+1 0,h - b b b • n n n 2 U n+1 0,h )v h dσ = ( U n k n , v h ) + 1 k n In (f, v h ) ds, ∀v h ∈ V i h , 3 k n (-U n+1 0,h + U n+1 1,h , w h ) + ã(U n+1 1,h , w h ) + Γ (Q n+1 0,h - b b b • n n n 2 U n+1 0,h )w h dσ = -( 3U n k n , w h ) + 6 k n In s -t n+1/2 k n (f, w h ) ds, ∀w h ∈ V i h , Γ Q n+1 0,h -b b b • n n nU n+1 0,h + p U n+1 0,h + q k n (U n+1 0,h + U n+1 1,h ) + q∇ Γ • (r r rU n+1 0,h ) ψ h dσ + Γ qs∇ Γ U n+1 0,h ∇ Γ ψ h dσ = Γ q U n k n ψ h dσ + 1 k n In Γ gψ h dσ ds, ∀ψ h ∈ W i,j h , Γ Q n+1 1,h -b b b • n n nU n+1 1,h + p U n+1 1,h + 3q k n (-U n+1 0,h + U n+1 1,h ) + q∇ Γ • (r r rU n+1 1,h ) ζ h dσ + Γ qs∇ Γ U n+1 1,h ∇ Γ ζ h dσ = - Γ q 3U n k n ζ h dσ + 6 k n In s -t n+1/2 k n Γ gζ h dσ ds, ∀ζ h ∈ W i,j h . Transmission terms. Let (U 0 i,h , Q 0 i,h ) ∈ P d (V i h , T i ) × P d ( W i,j h , T i ) be a given initial guess, for 1 ≤ i ≤ I. Then, at iteration k ≥ 1, we solve the subdomain problem in Ω i : d dt (I i U k i,h , v i,h ) i + ãi (U k i,h , v i,h ) i + Γi,j (Q k i,h - b b b i • n n n i 2 U k i,h )v i,h dσ = (P i f, v i,h ) i , in (0, T ), ∀v i,h ∈ V i h , Γi,j Q k i,h -b b b i • n n n i U k i,h + p i,j U k i,h ψ i,j h dσ + Γi,j q i,j (∂ t (I i U k i,h ) + ∇ Γi,j • (r r r i,j U k i,h ))ψ i,j h + q i,j s i,j ∇ Γi,j U k i,h ∇ Γi,j ψ i,j h dσ = P i gh ((U k-1 j,h , Q k-1 j,h ), ψ i,j h ), on (0, T ), ∀ψ i,j h ∈ W i,j h , j ∈ N i , (4.17)
with, for k ≥ 1,

gh ((U, Q), ψ) := Γi,j -Q + b b b j • n n n j U + p i,j U ψ dσ + Γi,j q i,j (∂ t (I j U ) + ∇ Γi,j • (r r r i,j U ))ψ + q i,j s i,j ∇ Γi,j U ∇ Γi,j ψ dσ. (4.18)
For the computation of the right-hand side in (4.17), we follow the same steps as in section 4.1.2, where we replace gk j,i with

g g g k j,i,h = gh ((U k-1 j,h , Q k-1 j,h ), ψ i,j 1,h ), ..., gh ((U k-1 j,h , Q k-1 j,h ), ψ i,j m i,j ,h ) t ,
and we replace Gi,k

m,α with Gi,k m,α,h ∈ W i,j h solution of d α=0 Gj,k m,α,h I j m ϕ j m,α ϕ j m,β ds = I j m g g g k j,i,h ϕ j m,β , β ∈ {0, ..., d}. (4.19)
For the computation of g g g k j,i,h , we write (U k-1 j,h ) |Γ i,j = n j,i l=1 u h j,l χ j,i l,h , and Q k-1 j,h = m j,i ℓ=1 z h j,ℓ ψ j,i ℓ,h , and introduce Q Q Q k-1 j,h = (z h j,1 , ..., z h j,m j,i ) t , U U U k-1 j,h = (u h j,1 , ..., u h j,n j,i ) t , and the projection matrices, for 1 ≤ k ≤ m i,j , 1 ≤ l ≤ m j,i , and 1 ≤ ℓ ≤ n j,i ,

( M i,j h ) k,l = Γi,j ψ i,j k,h ψ j,i l,h dσ, (M i,j h ) k,ℓ = Γi,j ψ i,j k,h χ j,i ℓ,h dσ, (M i,j b,h ) k,ℓ = Γi,j (-b b b j • n n n i + p i,j )ψ i,j k,h χ j,i ℓ,h dσ, (B i,j r,h ) k,ℓ = Γi,j ∇ Γi,j •(r r r i,j χ j,i ℓ,h )ψ i,j k,h dσ, (K i,j s,h ) k,ℓ = Γi,j q i,j s i,j ∇ Γi,j χ j,i ℓ,h ∇ Γi,j ψ i,j k,h dσ. Then g g g k j,i,h = -M i,j h Q Q Q k-1 j,h + ∂ t (I j M i,j h U U U k-1 j,h ) + (M i,j b,h + B i,j r,h + K i,j s,h )U U U k-1 j,h . (4.20)
The projection matrices are computed using the projection algorithm in [START_REF] Gander | A New Cement to Glue Nonconforming Grids with Robin Interface Conditions : The Finite Element Case[END_REF].

5. Numerical Results. We have implemented the algorithm with d = 1 and P 1 finite elements in space in each subdomain. Time windows are used in order to reduce the number of iterations of the algorithm. For the free parameters defining S i,j and Si,j , we chose r i,j to be the tangential component of the advection, s i,j the value of the diffusion in the domain Ω j . The optimized parameters p i,j and q i,j are constant along the interface. They correspond to a mean value of the parameters obtained by a numerical optimization of the convergence factor [START_REF] Gander | Schwarz waveform relaxation method for advectiondiffusion-reaction problems with discontinuous coefficients and non-matching grids[END_REF]. We first give an example of a multidomain solution with discontinuous variable diffusion, for two subdomains and one time window. The advection velocity is also discontinuous, taken normal to the interface in one subdomain, and tangential to the interface in the other subdomain. The latter case of a flow tangential to the interface is difficult when the interface conditions are not related to the convergence factor of the domain decomposition method (see for example [START_REF] Japhet | Méthode de décomposition de domaine et conditions aux limites artificielles en mécanique des fluides : méthode Optimisée d'Ordre 2 (OO2)[END_REF]).

The physical domain is Ω = (0, 1)×(0, 2), the final time is T = 1. The initial value is u 0 = 0.25e -100((x-0.55) 2 +(y-1.7) 2 ) and the right-hand side is f = 0. The domain Ω × (0, 2) is split into two subdomains Ω 1 = (0, 0.5) × (0, 2) and Ω 2 = (0. We analyze now the precision in time. The space mesh is conforming and the converged solution is such that the residual is smaller than 10 -8 . We compute a variational reference solution on a time grid with 4096 time steps. The nonconforming solutions are interpolated on the previous grid to compute the error. We start with a time grid with 128 time steps for the left domain and 94 time steps for the right domain. Thereafter the time step is divided by 2 several times. Figure 5.2 shows the norms of the error in L ∞ (I; L 2 (Ω i )) versus the number of refinements, for both subdomains. First we observe the order 2 in time for the nonconforming case. This fits the theoretical estimates, even though we have theoretical results only for Robin transmission conditions. Moreover, the error obtained in the nonconforming case, in the subdomain where the grid is finer, is nearly the same as the error obtained in the conforming finer case.

. Error between variational and DG-OSWR solutions versus the refinement in time

The computations are done using Order 2 transmission conditions. Indeed, the error between the multidomain and the variational solutions decrease much faster with the Order 2 transmissions conditions than with the Robin transmissions conditions as we can see in Figure 5.3, in the conforming case. We now consider the advection-diffusion equation with discontinuous porosity ω:

ω∂ t u + ∇ • (b b bu -ν∇u) = 0.
The physical domain is Ω = (0, 1) × (0, 2), the final time is T = 1.5. Ω is split into two subdomains. The interface Γ is parametrized with a Hermite polynomial ( 12 +((2s-1) 3 +2(2s-1) 2 +(2s-1))1 1 s≤ 1 2 +((2s-1) 3 -2(2s-1) 2 +(2s-1))1 1 s≥ 1 2 ), 0 < s < 1, see Figure 5 We first consider a conforming grid in space. The time step in Ω 1 is k 1 = 1/180, while in Ω 2 , k 2 = 1/100. In Figure 5.5, we observe, at final time T = 1.5, the approximate solution computed using ten time windows and 5 iterations in each time window. It is close to the variational solution computed in one time window on the conforming finer space-time grid as shown on the error, in Figure 5.6.

Conclusion.

We have proposed a new numerical method to solve parabolic equations with discontinuous coefficients. It relies on the splitting of the time interval into time windows, in which a few iterations of an optimized Schwarz waveform relaxation algorithm are performed by a discontinuous Galerkin method in time, with non conforming projection between space-time grids on the interfaces. We have shown theoretically in the Robin case that the method preserves the order of the discontinuous Galerkin method. Numerical estimates of the L 2 (I; L 2 (Ω i )) error and the H 1 error at final time have shown that the method preserves the order of the space nonconforming scheme as well. The analysis of the fully discrete scheme will be done in a further work.

Fig. 1 . 1 .

 11 Fig. 1.1. Left: decomposition with possible corners (Robin transmission conditions), right: decomposition in bands (Ventcell transmission conditions)

2. 1 .

 1 The local problem. The optimized Schwarz waveform relaxation algorithm relies on the resolution of the following initial boundary value problem in a domain O with boundary Γ:

  5, 1) × (0, 2). The reaction c is zero, the advection and diffusion coefficients are b b b 1 = (0, -1), ν 1 = 0.001 √ y, and b b b 2 = (-0.1, 0), ν 2 = 0.1 sin(xy). The mesh size over the interface and time step in Ω 1 are h 1 = 1/32 and k 1 = 1/128, while in Ω 2 , h 2 = 1/24 and k 2 = 1/94. In Figure 5.1, we observe, at final time T = 1, that the approximate solution computed using 3 iterations (right figure) is close to the variational solution computed in one time window on a time conforming finergrid (left figure).

Fig. 5 . 1 .

 51 Fig. 5.1. Variational (left) and nonconforming (right) DG-OSWR solutions

  Fig. 5.2. Error between variational and DG-OSWR solutions versus the refinement in time

Fig. 5 . 3 .

 53 Fig. 5.3. Convergence history for different transmission conditions

  .

4 .Fig. 5 . 4 .

 454 Fig. 5.4. Domains Ω 1 (left) and Ω 2 (right)

Fig. 5 . 5 .Fig. 5 . 6 .Fig. 5 . 8 . 2 Fig. 5 . 9 .Fig. 5 .

 5556582595 Fig. 5.5. DG-OSWR solution after 10 time windows and 5 iterations per window

The first two authors are partially supported by french ANR (COMMA)

are interpolated on the previous grid to compute the error. We start with a time grid with 120 time steps for the left domain and 26 time steps for the right domain and divide by 2 the time steps several times. Figure 5.7 shows the norms of the error in L ∞ (I; L 2 (Ω i )) versus the time steps, for both subdomains.

. Error curves versus the refinement in time

We observe the order 2 in time for the nonconforming case that fits the theoretical estimates. In Figure 5.7 we show also the norms of the error in L 2 (Ω i ) at final time t = T versus the time steps, for both subdomains. We observe the order 3 for the time nonconforming case. This corresponds to the superconvergence behavior described in [START_REF] Johnson | Adaptive finite element methods for parabolic problems ii: Optimal error estimates in l∞l 2 and l∞l∞[END_REF].

We now consider nonconforming grids in space as well. The mesh size and time step in Ω 1 are h 1 = 0.032 and k 1 = 1/120, while in Ω 2 , h 2 = 0.048 and k 2 = 1/26. In Figure 5.8 we observe, at final time T = 1.5, that the approximate solution computed using 5 iterations in one time window is close to the solution computed with the conformal in space grid in Figure 5.5, left. In Figure 5.9 and 5.10 we observe the precision versus the mesh size. The converged solution is such that the residual is smaller than 10 -12 . A variational reference solution is computed on a time grid with 960 time steps and a space grid with mesh size h = 3.5 10 -3 . The space-time nonconforming solutions are interpolated on the previous grid to compute the error. We start with a time grid with 60 time steps and a mesh size h 1 = 0.056 for the left domain and 20 time steps a mesh size h 2 = 0.11 for the right domain and divide by 2 the time step and mesh size several times. Figure 5.9 shows the norms of the error in L 2 (I; L 2 (Ω i )) versus the time steps, for both subdomains. We observe the order 2 for the nonconforming space-time case, even though we have theoretical results only for the time semi-discrete case. Figure 5.10 displays the norms of the error in L 2 (Ω i ) and in H 1 (Ω i ) at final time t = T versus the mesh size, for both subdomains. We observe the order 2 for the L 2 error, and the order 1 for the H 1 error for the nonconforming space-time case.