# Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 2: the variable coefficients case.

1. Introduction. The overall context of the present work is the
coupling between oceanic and atmospheric numerical models, in particular
for representing processes in which the interactions between both media
are of prime importance. The algorithms generally used to couple this
type of numerical models are often not fully correct from a mathematical
point of view. Indeed, they do not ensure a perfect consistency of the
fluxes exchanged at the air-sea interface [8]. In this context, the
long-term objective of our work is to derive alternative numerical
techniques ensuring such a consistency as well as to study their
possible impact on the physical results of coupled models.
Global-in-time optimized Schwarz methods (also called Schwarz waveform
relaxation methods) [4, 5] based on the concept of absorbing boundary
conditions [3] are particularly well suited for such problems. The
present study aims at finding efficient transmission conditions in the
case of the coupling between two diffusion equations representing the
turbulent vertical mixing in the planetary boundary layers near the
air-sea interface (see Section 3.1 for further details on the notion of
turbulent vertical mixing).

In the first part of this paper [9], we analytically derive optimized conditions in the case of a diffusion coefficient being constant in each medium but with a discontinuity through the interface. However, this provides only a simplified view of the true physics. The ocean and the atmosphere interact through various multi-scale physical processes that are usually hardly explicitly resolved by the spatio-temporal discretization. Because it is essential to account for the effect of the subgrid turbulent boundary layers on the resolved part of the flow, parameterization schemes were designed [7, 14]. Those schemes usually take the form of a turbulent mixing term with a spatially variable diffusion coefficient to account for local effects. Indeed, a parameterization with a constant diffusion originally introduced in [2] is now known to be naive. In this second part of the paper, we study the impact of this variability of the diffusion coefficients, in particular in the vicinity of the interface, on the convergence properties of the Schwarz algorithm. To our knowledge, the spatial variability of the coefficients has never been considered in the framework of Schwarz-like methods except in [10], where absorbing conditions are given for a one-dimensional stationary diffusion problem.

This paper is organized as follows. In the rest of this section, we briefly recall some theoretical aspects of optimized Schwarz methods necessary to follow the reasoning of the present study. In Section 2, we introduce a general methodology to analytically assess the impact of the spatial variability of the diffusion coefficients on the convergence of the Schwarz method. This method is applied first to a simple Dirichlet-Neumann algorithm and then to a more general Robin-Robin algorithm. Finally in Section 3, we illustrate the relevance of our approach by numerical results.

1.1. Model problem and Schwarz algorithm. The present study focuses on the coupling between two one-dimensional diffusion equations with variable coefficients. Consider two subdomains [[OMEGA].sub.1] =] - [L.sub.1], 0[ and [[OMEGA].sub.2] =]0, [L.sub.2][ with a common interface [GAMMA] = {x = 0}. The coupling problem reads

(1.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] correspond to the boundary conditions on the computational domain [OMEGA], and [F.sub.j] and [G.sub.j] are operators defining the interface conditions. Those operators must be designed to ensure a given consistency of the solution through [GAMMA]. In our study we require the equality of the subproblems solutions and of their normal fluxes.

In order to solve the coupling problem (1.1), we propose to implement a Schwarz algorithm with Robin-Robin interface conditions:

(1.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where k = 1, 2, ... is the iteration number and the initial guess [u.sup.0.sub.2] (0,t) is given. [[LAMBDA].sub.1] and [[LAMBDA].sub.1] are operators to be determined. As mentioned in [10], those operators can be either local or nonlocal.

1.2. Reminder of the framework in the case of constant (but discontinuous) diffusion coefficients. We briefly recall here some known results useful for the present study and detailed in [9]. The convergence study of the algorithm (1.2) with constant coefficients is performed by introducing the errors [e.sup.k.sub.j] = [u.sup.k.sub.j] - u* between the k-th iterate and the exact solution u* of the coupled problem. Using a Fourier transform in time (denoted for any function g [member of] [L.sup.2](R) by [??] := Fg), the partial differential equation [L.sub.j][e.sub.j] = 0 becomes an ordinary

differential equation [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is spatially constant here), whose characteristic roots are

[[sigma].sup.+.sub.j] = [square root of ([i[omega]/[D.sub.j]])], [[sigma].sup.-.sub.j] = - [[sigma].sup.+.sub.j] = - [square root of ([i[omega]/[D.sub.j]])].

It is then usually assumed that [L.sub.j] [right arrow] to and that [e.sub.j] tends to zero for x [right arrow] [infinity], which leads to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [alpha]([omega]) and [beta]([omega]) are determined to satisfy the boundary conditions. Finally, the convergence factor [rho] corresponding to the ratio between the errors at two successive iterations can be determined as a function of [[sigma].sup.[+ or -].sub.j], [D.sub.j], and [[lambda].sub.j] (the Fourier symbols of the operators [[LAMBDA].sub.j]):

(1.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We remark that in Fourier space, the following symbols

[[lambda].sup.opt.sub.1] = - [D.sub.2][[sigma].sup.-.sub.2] and [[lambda].sup.opt.sub.2] = [D.sub.1][[sigma].sup.+.sub.1]

lead to [rho] = 0, i.e., ensure convergence in two iterations However, the corresponding operators, which are called absorbing conditions, are nonlocal in time and therefore cannot be used in practical applications We thus need to look for a local approximation of these optimal operators. It was first suggested in [11] to use a low frequency approximation of the symbols based on a Taylor expansion around [omega] = 0 This results in effective transmission conditions only for [omega] being small. To obtain a more general approximation, efficient also for high frequencies, the so called optimized Schwarz methods (OSM) were introduced. The simplest version consists of approximating [[lambda].sup.opt.sub.1] and [[lambda].sup.opt.sub.2] by two constant values [[lambda].sup.0.sub.1] and [[lambda].sup.0.sub.2]: this corresponds to Robin interface conditions (also called zeroth order two-sided transmission conditions). The values for [[lambda].sup.0.sub.1] and [[lambda].sup.0.sub.2] are then determined by solving the optimization problem

(1.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In [9], this optimization problem is solved analytically for constant (and possibly discontinuous across [GAMMA]) diffusion coefficients. In this second part of our study, we complement the previous work [9] and discuss the effect of the spatial variability of the diffusion coefficients on the convergence speed and on the determination of the optimized conditions.

When the diffusion coefficient is spatially variable, the usual approach of determining the convergence factor is no longer straightforward. To circumvent this problem, we develop in the next section a methodology to analytically derive a convergence factor similar to (1.3) but including the spatial variability of the diffusion coefficients. Thanks to this new convergence factor, it will then be possible to find optimized values [[lambda].sup.0.sub.j] using (1.4). We expect a non-trivial effect of this variability on the convergence properties of the associated Schwarz algorithm. Indeed, in [10] it is shown for the stationary diffusion equation -[[partial derivative].sub.x](D(x)[[partial derivative].sub.x]u) = f that the absorbing conditions are given by Robin conditions with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. This result strongly suggests that it is not only the local values of the diffusion coefficient near the interface that have an impact on the parameters [[lambda].sub.j] but the whole profile D(x), x [member of] [OMEGA].

2. OSM for diffusion problems with spatially variable coefficients. As mentioned earlier, the diffusion coefficient may be spatially variable to account for local effects (e.g., in the turbulent boundary layers) within subdomains. In practical applications (like in oceanography or meteorology), diffusion coefficients are likely to vary by several orders of magnitude in the vertical direction (this point is further discussed in Section 3.1). This is the primary motivation to look for a methodology to analytically determine the convergence factor for non-constant diffusion coefficients defined on two non-overlapping subdomains. Throughout this study, we make the assumption that the diffusion profile does not vary with time.

2.1. Analytical determination of the shape of the errors. The first part of this section does not require any distinction between subdomains, so the j -subscripts are temporarily dropped. We denote by g(t) the function containing the information given by the neighboring subdomain, hence the problem under investigation is

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [lambda] being the Robin parameter we wish to determine to optimize the convergence speed. A Dirichlet condition is imposed at x = L, which corresponds in having [B.sub.1] = [B.sub.2] = I in (1.2) with I the identity map.

First, we notice that the method based on a Fourier analysis, commonly used to analytically determine the convergence factor, is less convenient for our model problem with variable coefficients. Indeed, in Fourier space we would obtain the ODE i[omega][??] - [[partial derivative].sub.x](D(x) [[partial derivative].sub.x][??]) = 0 for [??]. The study of this ODE appears to be at least as complicated as the original problem in physical space. This is why we propose to study directly the system (2.1). We transform this original problem with a homogeneous equation and nonhomogeneous boundary conditions into a problem with nonzero right-hand side but with homogeneous boundary conditions by searching for a solution of the form e(x, t) = [phi](x, t) + U(x, t) with [phi] being a lifting function satisfying the boundary conditions. The transformed problem reads

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The choice of [phi] is not unique. We choose this function as the solution of the problem (2.1) with a constant diffusion coefficient whose value is the value of D at x = 0, i.e., [phi] is the solution of

(2.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We then search for U(x,t) using a separation of variables U(x,t) = [summation over (n)][[PHI].sub.n](x)[T.sub.n](t). A

substitution in (2.2) leads to

[summation over (n)][T'.sub.n](t)[[PHI].sub.n](x) -[summation over (n)][T.sub.n](t)[[partial derivative].sub.x](D(x) [[partial derivative].sub.x] [[PHI].sub.n] (x)) = f (x,t),

where the right hand side is also expanded with respect to the functions j n(x) ,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The next step is to properly choose the functions [[PHI].sub.n]. An adequate choice would enable us to transform the PDE into ODEs for the unknown functions [[PHI].sub.n](x) and [T.sub.n](t). The natural choice is therefore to look for [[PHI].sub.n](x) as a solution of the following regular Sturm-Liouville (SL) problem

(2.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [c.sub.n] the eigenvalues of the SL operator. Such a choice leads to a family of functions [[PHI].sub.n](x) which are orthonormal for the scalar product (u, v) = [[integral].sup.L.sub.0] u(x)v(x)dx. The properties of regular SL problems are fully described in [1] or [6]. After some simple algebra, we find that a general solution of problem (2.1) is given by

(2.5) e(x, t) = f(x, t) + U(x, t) ,

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

satisfies (2.4), and [f.sub.n](t) satisfies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] with [??](x) = D(x) - D(0).

By formulating the solution of our problem using [??](x), we can properly separate the error into two parts corresponding to two different contributions: [phi](x, t) corresponds to the error for a constant coefficient D(0), and U(x, t) represents the error coming from the perturbations around D(0), namely [??](x).

We must now determine explicitly the function [??]. A straightforward way consists of using the continuous Fourier transform in time. By introducing the function [E.sub.w](x) = e[square root of ([i[omega]/D(0)]x] and by taking into account the boundary conditions at x = 0 and x = L, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is now possible to express the error (2.5) in the Fourier space. The functions [f.sub.n] are extended to zero for t < 0 and by the convolution theorem we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where H(t) is the Heaviside unit step function. The general form for [??](x,[omega]) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In practice it is usually assumed that the subdomains are unbounded (L [right arrow][infinity] to) to simplify the expression of the convergence factor and thus to simplify the optimization problem (1.4). Using this assumption, [??] becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which implies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

As a result of our study we get an expression for the error function in Fourier space that takes into account the spatial variability of the diffusion coefficient:

(2.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This error has been constructed for positive values of x, which can be identified as subdomain [[OMEGA].sub.2] following the notations introduced in Section. For negative x (i.e., on [[OMEGA].sub.1]), we obtain a very similar form:

(2.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the function h is the analog of the function g previously introduced.

The form of the error (2.6) suggests that the impact of the spatial variability of the diffusion coefficients will be primarily seen for low temporal frequencies. Indeed, the term [??](x) arising from the variability of the coefficient is weighted by [absolute value of [E.sub.[omega]](-x)]=[e.sup.-] [square root of ([omega]/2D(0)x)], making the effect of the variability negligible for large values of w but potentially significant for low frequencies. Moreover, we can draw the same conclusion for the variations with respect to x: when x is small (near the interface), [??](x) is weighted by a non-negligible number, while for x being large enough, [E.sub.[omega]] (-x) is very small.

2.2. Convergence factor of the Dirichlet-Neumann algorithm with spatially variable coefficients. So far we have established a general form of the errors propagating in each subdomain. We are now able to propose a formulation of the convergence speed for the global-intime Schwarz algorithm with spatially variable coefficients. Before dealing with the general Robin-Robin case, we intend to determine the convergence speed in a simpler Dirichlet-Neumann case, i.e., using the notations introduced in (1.1) for [G.sub.j] = I and [F.sub.j] = [D.sub.j](0)[[partial derivative]/ [[partial derivative].sub.x]] .

Moreover for the sake of practical convenience, we also try to find the expression of an "effective" value [D.sup.eff.sub.j] corresponding to a constant value which would have the same effect on the convergence speed as the non-constant diffusion profile [D.sub.j](x).

THEOREM 2.1 (Convergence factor with Dirichlet-Neumann transmission conditions). The convergence factor [[rho].sup.var.sub.DN] of the Schwarz algorithm (1.2) with Dirichlet-Neumann transmission conditions is

(2.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

(2.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and the eigenfunctions [[PHI].sub.n,j] and eigenvalues [c.sub.n,j] being solutions of the Sturm-Liouville problem (2.4).

Proof. Hereafter we use again the subscripts j to characterize both subdomains, and we use the function [E.sub.[omega],j] (x) defined above that plays the same role as the function [E.sub.[omega]] previously defined. A derivation very similar to what has been done in Section 2.1 but with a Dirichlet boundary condition instead of a Robin boundary condition leads to

(2.10) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [??]([omega]) = [[??].sub.1](0, [omega]) and where the functions [[PHI].sub.n,2] are defined by a SL problem similar to (2.4) but again with a Dirichlet condition instead of a Robin condition. On [[OMEGA].sub.1], we have (by simply taking [lambda] = 0 in the derivation of Section 2.1):

(2.11) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and where the functions [[PHI].sub.n,1] are defined by a SL problem similar to (2.4) with a homogeneous Neumann condition at x = 0. The multiplicative Schwarz algorithm with Dirichlet-Neumann conditions is obtained by replacing [[??].sub.2] (respectively [??]) by [[??].sup.k.sub.2], (respectively [[??].sup.k.sub.2] (0, [omega])) in (2.10), and [[??].sub.1] (respectively [??]) by [[??].sup.k.sub.2] (respectively [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in (2.11). Therefore we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then, if we define a convergence factor by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

the previous relations lead to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[rho].sup.cst.sub.DN] = [quare root of ([D.sub.2](0)/[D.sub.1](0) is the convergence factor obtained in the case of constant diffusion coefficients (see [12]) and [[??].sub.DN] is given in (2.9). ?

Theorem 2.1 shows that the convergence factor [[rho].var.cst.sub.DN] naturally appears as the product of the convergence factor with constant coefficients (the surface values) and a term coming from the spatial variability of the diffusion coefficient on [[OMEGA].sub.1] and [[OMEGA].sub.1].

Starting from Equation (2.8), we can suggest two "effective" constant values for [D.sub.1] and [D.sub.2]. These (spatially constant) values have a similar effect on the convergence speed as the non-constant vertical profiles [D.sub.1](x) and [D.sub.2](x). They satisfy [[rho].sup.cst.sub.DN] = [square root of ([D.sup.eff.sub.2]([omega])/ [D.sup.eff.sub.2]([omega])] with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is worth mentioning that, due to the variability of the coefficients, the convergence factor is a function of the time frequency [omega], whereas this dependency does not exist with constant coefficients. Some examples of convergence factors [[rho].var.cst.sub.DN] are given in Section 3.2. Note that in the case [omega] [right arrow] 0, we get [D.sup.eff.sub.1] [right arrow] [D.sub.1](0), while

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The effect of the variability of the coefficient in the subdomain with a Neumann condition asymptotically vanishes. This is, however, not the case for the subdomain [[OMEGA].sub.2] with Dirichlet conditions ([D.sup.eff.sub.2] ([omega] [right arrow] 0) [not equal to] [D.sub.2](0)). This result shows that, when a Dirichlet-Neumann algorithm is used, [[rho].sup.cst.sub.DN] < 1 does not necessarily imply that [[rho].sup.var.sub.DN] < 1. In other words, the fact that the algorithm with constant coefficients (the interface values) theoretically converges does not ensure that the algorithm with variable coefficients will. Indeed,

(2.12)[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

2.3. Convergence factor of the Robin-Robin algorithm with spatially variable coefficients. In this section we determine the convergence factor [[rho].sup.var.sub.RR] for the more general case of Robin-Robin interface conditions.

THEOREM 2.2 (Convergence factor with Robin-Robin transmission conditions). The convergence factor [[rho].sup.var.sub.RR] of the Schwarz algorithm (1.2) with Robin-Robin transmission conditions is

(2.13) [[rho].sup.var.sub.RR] = [absolute value of [([[lambda].sub.1] + [[lambda].sub.2])[K.sub.1] - 1] [([[lambda].sub.1] + [[lambda].sub.2])[K.sub.2] - 1]],

with [[lambda].sub.j] being the Fourier symbol of the operator [[LAMBDA].sub.j] in (1.2) and

(2.14) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and the eigenfunctions [[PHI].sub.n,j] and eigenvalues [c.sub.n,j] are solutions of the Sturm-Liouville problem (2.4).

Proof. Thanks to (2.6) and (2.7), we can express [[??].sub.1] and [[??].sub.2] in a compact form for the iterate k as

(2.15) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [K.sub.j] is given in (2.14). The problem on the interface x = 0 is given by the relations

(2.16) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and by combining (2.15) and (2.16), we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

By substituting those expressions in (2.16), we finally get a relation linking [??] and [??]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which leads to an expression for the convergence factor

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We note that this expression of the convergence factor is consistent with the expression (1.3) obtained in the case of constant (but discontinuous) coefficients. Indeed, if we set [[??].sub.1](x) = [[??].sub.2] (x) =0 in (2.13), we then have [K.sub.j] = 1/[square root of (i[omega][D.sub.j](0), which leads to (1.3) because [D.sub.j][[sigma].sup.[+ or -].sub.j] = [+ or -][square root of i[omega][D.sub.j](0)]. A convenient formulation of [[rho].sup.var.sub.RR], without complex numbers can be found in Appendix A. To conclude this section, we look at the asymptotic beha ior of [[rho].sup.var.sub.RR] and we can easily find that

[[rho].sup.var.sub.RR] ([omega] [right arrow] 0) = [[rho].sup.var.sub.RR] ([omega] [right arrow] [infinity]) [right arrow] 1,

which shows that the effect of the variability of the diffusion coefficients asymptotically vanishes when a Robin-Robin algorithm is used.

3. Numerical results. In this section we verify numerically the validity of the theoretical results presented in Section 2. To do this, we first briefly describe the rationale for the spatial variability of the diffusion coefficient and provide a typical profile which we will use for the numerical tests. Then we design a few experiments to illustrate the relevance of our theoretical results.

3.1. Planetary boundary layer turbulence. Unlike boundary layers in many engineering flows, the atmospheric and oceanic planetary boundary layers are almost always turbulent and cannot be explicitly resolved due to the insufficient vertical resolution in computational models. The numerical representation of those layers thus relies on the Reynolds decomposition: the flow is split into a mean (resolved) part <u> and a fluctuating (subgrid) part u' (where u can either represent a velocity component or an active tracer). When this decomposition is applied to nonlinear (advective) terms, this gives rise to additional terms and hence to a closure problem. The dominant expression in the turbulent boundary layers arising from the Reynolds decomposition is the divergence of the vertical part of (u 'w ') (where w denotes the vertical component of the velocity). Typically, this turbulent vertical flux is expressed as a function of the mean (resolved) part of the flow by using the down-gradient assumption, (u 'w ') = -D(x)[[partial derivative].sub.x] (u), where D(x) is the so-called eddy diffusivity or eddy-viscosity if u represents a velocity. This assumption explains why a one-dimensional diffusion equation, like the one studied in the present paper, is generally sufficient to locally represent the turbulent mixing in the boundary layers. The eddy diffusivity D(x) is defined to allow the flow to make the transition between its surface (the air-sea interface) and its interior (below the boundary layer) properties. This is the reason why D(x) exhibits a strong spatial variability. In this context, several ways to specify the coefficient D(x) have been proposed. The formulation most commonly used in the state-of-the-art numerical models can be found in [7] and [14]. Those formulations define the eddy diffusivity as

(3.1)[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [h.sub.bl] the thickness of the boundary layer (depending on the state of the flow) and A a parameter setting the intensity of the mixing (note that D(x) is continuous and differentiable at x = [h.sub.bl]). Throughout this section we assume that D(x) is given by (3.1), and a typical profile for A = 0.5 m/s and [h.sub.bl] = 150 m is given in Figure 3.1.

[FIGURE 3.1 OMITTED]

In the remainder of this section we study first a Dirichlet-Neumann algorithm and then a Robin-Robin algorithm. We define spatially variable coefficients with AT =0.1 m/s (respectively [A.sub.2] = 0.5 m/s) and [h.sub.bl,1] = 50 m (respectively [h.sub.bl,2] = 150 m) on [[OMEGA].sub.1] (respectively [[OMEGA].sub.2]). The values of [v.sub.1] and [v.sub.2] (corresponding to the surface values [D.sub.1](0) and [D.sub.2] (0)) are chosen to be the same as the values used in [9] in the constant coefficient case. If we introduce [gamma] = [square root of ([[v.sub.2]/[v.sub.1]])], we investigate the two cases [gamma] =10 and [gamma] =[square root of ([square root of (10)]] with [v.sub.2] = 0.5 [m.sup.2]/s (the value of [v.sup.1] is adjusted depending on the value of [gamma]). Those various parameter values lead to diffusion profiles that can be found in the atmospheric and oceanic boundary layers. The discretization of the problem, the computational grid, as well as the initial conditions are described in [9, Section 5]. We use [DELTA]t = 100 s and a random initial guess on the interface so that it contains a wide range of the temporal frequencies that can be resolved by the computational grid.

3.2. Test case 1: Dirichlet-Neumann. The analytical convergence factor [[rho].sup.var.sub.DN] ([omega]) in Equation (2.8) is shown for different values of [gamma] in Figure 3.2. The eigenvalues [c.sub.n] and eigenfunctions [[PHI].sub.n] are computed numerically on the same computational grid as the model problem. We remark that depending on the jump in the coefficients through the interface, the spatial variability of the diffusivities either tend to accelerate the convergence speed (for [gamma] = [square root of ([square root of (10)]]) or to slow it down (for [gamma] = 10) compared to the convergence speed obtained with constant coefficients. As expected, the convergence factor for spatially variable coefficients is no longer independent of [omega], and for low frequencies we get a significant departure from the convergence rate of the algorithm with constant coefficients. The trend seen in the convergence factor [[rho].sup.var.sub.DN] ([omega]) determined at a continuous level is confirmed by the numerical results (Figure 3.2, top). These results as well as the asymptotic expression (2.12) call for caution when we use a Dirichlet-Neumann algorithm with spatially variable coefficients because it can lead to significantly different performances compared to the one obtained with constant coefficients. We can expect a Robin-Robin type algorithm to provide a more robust alternative thanks to the tuning of the [[lambda].sub.j] parameters.

[FIGURE 3.2 OMITTED]

3.3. Test case 2: Robin-Robin. Here, we denote by [[lambda].sup.cst.sub.j] the optimal Robin parameters obtained using the analytical results found in [9] for constant coefficients. We consider that these constant coefficients are the interface values [D.sub.j](0). Moreover, we denote by [[lambda].sup.var.sub.j] the Robin parameters optimized by solving numerically the problem (1.4) with the convergence factor [[rho].sup.var.sub.RR] as given in (2.13). This optimization is done using the Rosenbrock method [13] and by taking the parameters [[lambda].sup.cst.sub.j] to initialize the algorithm. We see from Figure 3.3 that the use of the parameters [[lambda].sup.var.sub.j] provide slightly better convergence properties compared to the parameters [[lambda].sup.cst.sub.j], regardless of the value of [gamma]. As for the Dirichlet-Neumann algorithm, we can check that our analytical study at the continuous level provides a convergence factor [[rho].sup.var.sub.RR] representative of the behavior of the algorithm at a discrete level (Figure 3.3, bottom). From Figure 3.3 (bottom left), we also see that the way we initialize the algorithm (with a random initial guess for [u.sup.0.sub.2] (0, t), t [member of] [0,T]) leads to the generation of a large range of temporal frequencies and more particularly low frequencies slowing down the convergence speed of the simulation using the parameters [[lambda].sup.cst.sub.j], although the latter provide a faster convergence than the parameters [[lambda].sup.var.sub.j] for most of the frequency spectrum. For our model problem, the use of the parameters [[lambda].sup.var.sub.j] provides a relatively modest improvement over the [[lambda].sup.cst.sub.j] parameters. However, in general, this statement has to be mitigated because if we consider [h.sub.bl,2] = 10 m instead of [h.sub.bl,2] = 150 m, we see in Figure 3.4 that the parameters obtained through an optimization of [[rho].sup.var.sub.RR] are clearly superior to the parameters [[lambda].sup.cst.sub.j]. For this case, we show in Figure 3.5 the asymptotic behavior of the optimized convergence rate and the associated Robin parameters [[lambda].sup.var.sub.j]. We numerically show that thanks to the theoretical work presented earlier in the paper, we are able to bring the proper adjustments to the Robin parameters so that our algorithm asymptotically maintains a good efficiency even in the presence of spatially variable coefficients.

[FIGURE 3.3 OMITTED]

4. Conclusion. In this paper we present and analyze a new approach to study the convergence properties of a global-in-time Schwarz algorithm in the case of a one-dimensional diffusion problem with spatially variable diffusion coefficients. We analytically derive an expression for the evolution of the errors of such an algorithm with respect to the iterates. Thanks to our formulation, we are able to gain a better understanding of the behavior of the associated convergence factor. We exhibit some interesting features that were not shown by the usual convergence studies with constant diffusion coefficients. We put particular emphasis on the fact that for low temporal frequencies, it can be a very inaccurate assumption to replace a variable diffusion coefficient by its constant interface value in the convergence study. Moreover, we also show that depending on the type of algorithm under consideration (Dirichlet-Neumann or Robin-Robin) the variability of the coefficients may have more or less impact on the asymptotic convergence properties. To be more attractive for practical applications, our approach requires further developments by performing an accurate study of the eigenvalue problems to improve our knowledge of the behavior of these eigenvalues with respect to the perturbations of the diffusion profiles.

[FIGURE 3.4 OMITTED]

[FIGURE 3.5 OMITTED]

Appendix A. Determination of the convergence factor in the case of variable coefficients. We recall (2.13):

(A.1) [rho] = [absolute value of [([[lambda].sub.1] + [lambda].sub.2])[K.sub.1] - 1] [([lambda].sub.1] + [lambda].sub.2])[K.sub.2] - 1]],

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(A.1) can be rewritten as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In order to determine the real and imaginary parts of [K.sub.1] and [K.sub.2], we can decompose each term appearing in the preceding expressions:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thanks to these equalities, we can recast Kj into the following form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and by noting that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we obtain

[K.sub.1] = ([e.sub.1](1 - [g.sub.1]) + [f.sub.1][h.sub.1]) + i([f.sub.1](1 - [g.sub.1]) - [e.sub.1][h.sub.1]),

[K.sub.2] = ([e.sub.2](1 - [g.sub.2]) + [f.sub.1][h.sub.2]) + i([f.sub.2](1 - [g.sub.2]) - [e.sub.1][h.sub.2]).

Hence, thanks to (A.2), this yields an expression for the convergence factor p without complex numbers.

Acknowledgments. This research was partially supported by the ANR project "COMMA" and by the INRIA project-team MOISE. The authors are thankful for the comments and suggestions of two anonymous reviewers, which helped to improve the clarity of this manuscript.

REFERENCES

[1] P. B. BAILEY, W. N. EVERITT, AND A. ZETTL, Regular and singular Sturm-Liouville problems with coupled boundary conditions, Proc. Roy. Soc. Edinburgh Sect. A, 126 (1996), pp. 505-514.

[2] V. W. EKMAN, On the influence of'the earth's rotation on ocean currents, Ark. Mat. Astron. Fys., 2 (1905), pp. 1-53.

[3] B. ENGQUIST AND A. MAJDA, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp., 31 (1977), pp. 629-651.

[4] M.J. GANDER ANDL. HALPERN, Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems, SIAM J. Numer. Anal., 45 (2007), pp. 666-697.

[5] M. J. GANDER, L. HALPERN, AND F. NATAF, Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation, in Eleventh International Conference on Domain Decomposition Methods London 1998, C. H. Lai, P. E. Bjorstad, M. Cross, and O. Widlund, eds., DDM.org, Augsburg, 1999, pp. 27-36.

[6] Q. KONG AND A. ZETTL, Eigenvalues ofregular Sturm-Liouville problems, J. Differential Equations, 131 (1996), pp. 1-19.

[7] W. G. LARGE, J. C. MCWILLIAMS, AND S. C. DONEY, Oceanic vertical mixing: A review and a model with a nonlocal boundarylayer parameterization, Rev. Geophys., 32 (1994), pp. 363-403.

[8] F. LEMARIE, Algorithmes de Schwarz et couplage ocean-atmosphere, Ph.D. Thesis, Laboratoire Jean Kuntz mann, Joseph Fourier University, Grenoble, 2008.

[9] F. LEMARlE, L. DEBREU, AND E. BLAYO, Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 1: the constant coefficients case, Electron. Trans. Numer. Anal., 40 (2013), pp. 148-169. http://etna.math.kent.edu/vol.40.2 013/pp14 8-16 9.dir

[10] P.-L. LIONS, On the Schwarz alternating method. III. A variant for nonoverlapping subdomains, in Third International Symposium on Domain Decomposition Methods for Partial Differential Equations, T. Chan, R. Glowinski, J. Periaux, and O. Widlund, eds., SIAM, Philadelphia, 1990, pp. 202-223.

[11] F. NATAF, F. ROGIER, AND E. DE STURLER, Optimal interface conditions for domain decomposition methods, Internal Report 301, CMAP, Ecole Polytechnique, Palaiseau, 1994.

[12] A. QUARTERONI AND A. VALLI, Domain Decomposition Methods for Partial Differential Equations, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, 1999.

[13] H. H. ROSENBROCK, An automatic method for finding the greatest or least value of a function, Comput. J., 3 (1960), pp. 175-184.

[14] I. B. TROEN AND L. MAHRT, A simple model of the atmospheric boundary layer; sensitivity to surface evaporation, Boundary-Layer Meteorol., 37 (1986), pp. 129-148.

FLORIAN LEMARIE[dagger], LAURENT DEBREU[dagger], AND ERIC BLAYO[double dagger]

* Received September 29, 2010. Accepted March 11, 2013. Published online on July 5, 2013. Recommended by Martin Gander.

[dagger] INRIA Grenoble Rhone-Alpes, Montbonnot, 38334 Saint Ismier Cedex, France and Jean Kuntzmann Laboratory, BP 53, 38041 Grenoble Cedex 9, France, ({florian. lemarie, laurent. debreu}@inria .fr).

[double dagger] University of Grenoble, Jean Kuntzmann Laboratory, BP 53, 38041 Grenoble Cedex 9, France (eric.blayo@imag.fr).

In the first part of this paper [9], we analytically derive optimized conditions in the case of a diffusion coefficient being constant in each medium but with a discontinuity through the interface. However, this provides only a simplified view of the true physics. The ocean and the atmosphere interact through various multi-scale physical processes that are usually hardly explicitly resolved by the spatio-temporal discretization. Because it is essential to account for the effect of the subgrid turbulent boundary layers on the resolved part of the flow, parameterization schemes were designed [7, 14]. Those schemes usually take the form of a turbulent mixing term with a spatially variable diffusion coefficient to account for local effects. Indeed, a parameterization with a constant diffusion originally introduced in [2] is now known to be naive. In this second part of the paper, we study the impact of this variability of the diffusion coefficients, in particular in the vicinity of the interface, on the convergence properties of the Schwarz algorithm. To our knowledge, the spatial variability of the coefficients has never been considered in the framework of Schwarz-like methods except in [10], where absorbing conditions are given for a one-dimensional stationary diffusion problem.

This paper is organized as follows. In the rest of this section, we briefly recall some theoretical aspects of optimized Schwarz methods necessary to follow the reasoning of the present study. In Section 2, we introduce a general methodology to analytically assess the impact of the spatial variability of the diffusion coefficients on the convergence of the Schwarz method. This method is applied first to a simple Dirichlet-Neumann algorithm and then to a more general Robin-Robin algorithm. Finally in Section 3, we illustrate the relevance of our approach by numerical results.

1.1. Model problem and Schwarz algorithm. The present study focuses on the coupling between two one-dimensional diffusion equations with variable coefficients. Consider two subdomains [[OMEGA].sub.1] =] - [L.sub.1], 0[ and [[OMEGA].sub.2] =]0, [L.sub.2][ with a common interface [GAMMA] = {x = 0}. The coupling problem reads

(1.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] correspond to the boundary conditions on the computational domain [OMEGA], and [F.sub.j] and [G.sub.j] are operators defining the interface conditions. Those operators must be designed to ensure a given consistency of the solution through [GAMMA]. In our study we require the equality of the subproblems solutions and of their normal fluxes.

In order to solve the coupling problem (1.1), we propose to implement a Schwarz algorithm with Robin-Robin interface conditions:

(1.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where k = 1, 2, ... is the iteration number and the initial guess [u.sup.0.sub.2] (0,t) is given. [[LAMBDA].sub.1] and [[LAMBDA].sub.1] are operators to be determined. As mentioned in [10], those operators can be either local or nonlocal.

1.2. Reminder of the framework in the case of constant (but discontinuous) diffusion coefficients. We briefly recall here some known results useful for the present study and detailed in [9]. The convergence study of the algorithm (1.2) with constant coefficients is performed by introducing the errors [e.sup.k.sub.j] = [u.sup.k.sub.j] - u* between the k-th iterate and the exact solution u* of the coupled problem. Using a Fourier transform in time (denoted for any function g [member of] [L.sup.2](R) by [??] := Fg), the partial differential equation [L.sub.j][e.sub.j] = 0 becomes an ordinary

differential equation [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is spatially constant here), whose characteristic roots are

[[sigma].sup.+.sub.j] = [square root of ([i[omega]/[D.sub.j]])], [[sigma].sup.-.sub.j] = - [[sigma].sup.+.sub.j] = - [square root of ([i[omega]/[D.sub.j]])].

It is then usually assumed that [L.sub.j] [right arrow] to and that [e.sub.j] tends to zero for x [right arrow] [infinity], which leads to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [alpha]([omega]) and [beta]([omega]) are determined to satisfy the boundary conditions. Finally, the convergence factor [rho] corresponding to the ratio between the errors at two successive iterations can be determined as a function of [[sigma].sup.[+ or -].sub.j], [D.sub.j], and [[lambda].sub.j] (the Fourier symbols of the operators [[LAMBDA].sub.j]):

(1.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We remark that in Fourier space, the following symbols

[[lambda].sup.opt.sub.1] = - [D.sub.2][[sigma].sup.-.sub.2] and [[lambda].sup.opt.sub.2] = [D.sub.1][[sigma].sup.+.sub.1]

lead to [rho] = 0, i.e., ensure convergence in two iterations However, the corresponding operators, which are called absorbing conditions, are nonlocal in time and therefore cannot be used in practical applications We thus need to look for a local approximation of these optimal operators. It was first suggested in [11] to use a low frequency approximation of the symbols based on a Taylor expansion around [omega] = 0 This results in effective transmission conditions only for [omega] being small. To obtain a more general approximation, efficient also for high frequencies, the so called optimized Schwarz methods (OSM) were introduced. The simplest version consists of approximating [[lambda].sup.opt.sub.1] and [[lambda].sup.opt.sub.2] by two constant values [[lambda].sup.0.sub.1] and [[lambda].sup.0.sub.2]: this corresponds to Robin interface conditions (also called zeroth order two-sided transmission conditions). The values for [[lambda].sup.0.sub.1] and [[lambda].sup.0.sub.2] are then determined by solving the optimization problem

(1.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In [9], this optimization problem is solved analytically for constant (and possibly discontinuous across [GAMMA]) diffusion coefficients. In this second part of our study, we complement the previous work [9] and discuss the effect of the spatial variability of the diffusion coefficients on the convergence speed and on the determination of the optimized conditions.

When the diffusion coefficient is spatially variable, the usual approach of determining the convergence factor is no longer straightforward. To circumvent this problem, we develop in the next section a methodology to analytically derive a convergence factor similar to (1.3) but including the spatial variability of the diffusion coefficients. Thanks to this new convergence factor, it will then be possible to find optimized values [[lambda].sup.0.sub.j] using (1.4). We expect a non-trivial effect of this variability on the convergence properties of the associated Schwarz algorithm. Indeed, in [10] it is shown for the stationary diffusion equation -[[partial derivative].sub.x](D(x)[[partial derivative].sub.x]u) = f that the absorbing conditions are given by Robin conditions with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. This result strongly suggests that it is not only the local values of the diffusion coefficient near the interface that have an impact on the parameters [[lambda].sub.j] but the whole profile D(x), x [member of] [OMEGA].

2. OSM for diffusion problems with spatially variable coefficients. As mentioned earlier, the diffusion coefficient may be spatially variable to account for local effects (e.g., in the turbulent boundary layers) within subdomains. In practical applications (like in oceanography or meteorology), diffusion coefficients are likely to vary by several orders of magnitude in the vertical direction (this point is further discussed in Section 3.1). This is the primary motivation to look for a methodology to analytically determine the convergence factor for non-constant diffusion coefficients defined on two non-overlapping subdomains. Throughout this study, we make the assumption that the diffusion profile does not vary with time.

2.1. Analytical determination of the shape of the errors. The first part of this section does not require any distinction between subdomains, so the j -subscripts are temporarily dropped. We denote by g(t) the function containing the information given by the neighboring subdomain, hence the problem under investigation is

(2.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [lambda] being the Robin parameter we wish to determine to optimize the convergence speed. A Dirichlet condition is imposed at x = L, which corresponds in having [B.sub.1] = [B.sub.2] = I in (1.2) with I the identity map.

First, we notice that the method based on a Fourier analysis, commonly used to analytically determine the convergence factor, is less convenient for our model problem with variable coefficients. Indeed, in Fourier space we would obtain the ODE i[omega][??] - [[partial derivative].sub.x](D(x) [[partial derivative].sub.x][??]) = 0 for [??]. The study of this ODE appears to be at least as complicated as the original problem in physical space. This is why we propose to study directly the system (2.1). We transform this original problem with a homogeneous equation and nonhomogeneous boundary conditions into a problem with nonzero right-hand side but with homogeneous boundary conditions by searching for a solution of the form e(x, t) = [phi](x, t) + U(x, t) with [phi] being a lifting function satisfying the boundary conditions. The transformed problem reads

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The choice of [phi] is not unique. We choose this function as the solution of the problem (2.1) with a constant diffusion coefficient whose value is the value of D at x = 0, i.e., [phi] is the solution of

(2.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We then search for U(x,t) using a separation of variables U(x,t) = [summation over (n)][[PHI].sub.n](x)[T.sub.n](t). A

substitution in (2.2) leads to

[summation over (n)][T'.sub.n](t)[[PHI].sub.n](x) -[summation over (n)][T.sub.n](t)[[partial derivative].sub.x](D(x) [[partial derivative].sub.x] [[PHI].sub.n] (x)) = f (x,t),

where the right hand side is also expanded with respect to the functions j n(x) ,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The next step is to properly choose the functions [[PHI].sub.n]. An adequate choice would enable us to transform the PDE into ODEs for the unknown functions [[PHI].sub.n](x) and [T.sub.n](t). The natural choice is therefore to look for [[PHI].sub.n](x) as a solution of the following regular Sturm-Liouville (SL) problem

(2.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [c.sub.n] the eigenvalues of the SL operator. Such a choice leads to a family of functions [[PHI].sub.n](x) which are orthonormal for the scalar product (u, v) = [[integral].sup.L.sub.0] u(x)v(x)dx. The properties of regular SL problems are fully described in [1] or [6]. After some simple algebra, we find that a general solution of problem (2.1) is given by

(2.5) e(x, t) = f(x, t) + U(x, t) ,

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

satisfies (2.4), and [f.sub.n](t) satisfies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] with [??](x) = D(x) - D(0).

By formulating the solution of our problem using [??](x), we can properly separate the error into two parts corresponding to two different contributions: [phi](x, t) corresponds to the error for a constant coefficient D(0), and U(x, t) represents the error coming from the perturbations around D(0), namely [??](x).

We must now determine explicitly the function [??]. A straightforward way consists of using the continuous Fourier transform in time. By introducing the function [E.sub.w](x) = e[square root of ([i[omega]/D(0)]x] and by taking into account the boundary conditions at x = 0 and x = L, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is now possible to express the error (2.5) in the Fourier space. The functions [f.sub.n] are extended to zero for t < 0 and by the convolution theorem we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where H(t) is the Heaviside unit step function. The general form for [??](x,[omega]) is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In practice it is usually assumed that the subdomains are unbounded (L [right arrow][infinity] to) to simplify the expression of the convergence factor and thus to simplify the optimization problem (1.4). Using this assumption, [??] becomes

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which implies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

As a result of our study we get an expression for the error function in Fourier space that takes into account the spatial variability of the diffusion coefficient:

(2.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This error has been constructed for positive values of x, which can be identified as subdomain [[OMEGA].sub.2] following the notations introduced in Section. For negative x (i.e., on [[OMEGA].sub.1]), we obtain a very similar form:

(2.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where the function h is the analog of the function g previously introduced.

The form of the error (2.6) suggests that the impact of the spatial variability of the diffusion coefficients will be primarily seen for low temporal frequencies. Indeed, the term [??](x) arising from the variability of the coefficient is weighted by [absolute value of [E.sub.[omega]](-x)]=[e.sup.-] [square root of ([omega]/2D(0)x)], making the effect of the variability negligible for large values of w but potentially significant for low frequencies. Moreover, we can draw the same conclusion for the variations with respect to x: when x is small (near the interface), [??](x) is weighted by a non-negligible number, while for x being large enough, [E.sub.[omega]] (-x) is very small.

2.2. Convergence factor of the Dirichlet-Neumann algorithm with spatially variable coefficients. So far we have established a general form of the errors propagating in each subdomain. We are now able to propose a formulation of the convergence speed for the global-intime Schwarz algorithm with spatially variable coefficients. Before dealing with the general Robin-Robin case, we intend to determine the convergence speed in a simpler Dirichlet-Neumann case, i.e., using the notations introduced in (1.1) for [G.sub.j] = I and [F.sub.j] = [D.sub.j](0)[[partial derivative]/ [[partial derivative].sub.x]] .

Moreover for the sake of practical convenience, we also try to find the expression of an "effective" value [D.sup.eff.sub.j] corresponding to a constant value which would have the same effect on the convergence speed as the non-constant diffusion profile [D.sub.j](x).

THEOREM 2.1 (Convergence factor with Dirichlet-Neumann transmission conditions). The convergence factor [[rho].sup.var.sub.DN] of the Schwarz algorithm (1.2) with Dirichlet-Neumann transmission conditions is

(2.8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where

(2.9) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and the eigenfunctions [[PHI].sub.n,j] and eigenvalues [c.sub.n,j] being solutions of the Sturm-Liouville problem (2.4).

Proof. Hereafter we use again the subscripts j to characterize both subdomains, and we use the function [E.sub.[omega],j] (x) defined above that plays the same role as the function [E.sub.[omega]] previously defined. A derivation very similar to what has been done in Section 2.1 but with a Dirichlet boundary condition instead of a Robin boundary condition leads to

(2.10) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [??]([omega]) = [[??].sub.1](0, [omega]) and where the functions [[PHI].sub.n,2] are defined by a SL problem similar to (2.4) but again with a Dirichlet condition instead of a Robin condition. On [[OMEGA].sub.1], we have (by simply taking [lambda] = 0 in the derivation of Section 2.1):

(2.11) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and where the functions [[PHI].sub.n,1] are defined by a SL problem similar to (2.4) with a homogeneous Neumann condition at x = 0. The multiplicative Schwarz algorithm with Dirichlet-Neumann conditions is obtained by replacing [[??].sub.2] (respectively [??]) by [[??].sup.k.sub.2], (respectively [[??].sup.k.sub.2] (0, [omega])) in (2.10), and [[??].sub.1] (respectively [??]) by [[??].sup.k.sub.2] (respectively [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in (2.11). Therefore we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then, if we define a convergence factor by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

the previous relations lead to

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[rho].sup.cst.sub.DN] = [quare root of ([D.sub.2](0)/[D.sub.1](0) is the convergence factor obtained in the case of constant diffusion coefficients (see [12]) and [[??].sub.DN] is given in (2.9). ?

Theorem 2.1 shows that the convergence factor [[rho].var.cst.sub.DN] naturally appears as the product of the convergence factor with constant coefficients (the surface values) and a term coming from the spatial variability of the diffusion coefficient on [[OMEGA].sub.1] and [[OMEGA].sub.1].

Starting from Equation (2.8), we can suggest two "effective" constant values for [D.sub.1] and [D.sub.2]. These (spatially constant) values have a similar effect on the convergence speed as the non-constant vertical profiles [D.sub.1](x) and [D.sub.2](x). They satisfy [[rho].sup.cst.sub.DN] = [square root of ([D.sup.eff.sub.2]([omega])/ [D.sup.eff.sub.2]([omega])] with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is worth mentioning that, due to the variability of the coefficients, the convergence factor is a function of the time frequency [omega], whereas this dependency does not exist with constant coefficients. Some examples of convergence factors [[rho].var.cst.sub.DN] are given in Section 3.2. Note that in the case [omega] [right arrow] 0, we get [D.sup.eff.sub.1] [right arrow] [D.sub.1](0), while

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The effect of the variability of the coefficient in the subdomain with a Neumann condition asymptotically vanishes. This is, however, not the case for the subdomain [[OMEGA].sub.2] with Dirichlet conditions ([D.sup.eff.sub.2] ([omega] [right arrow] 0) [not equal to] [D.sub.2](0)). This result shows that, when a Dirichlet-Neumann algorithm is used, [[rho].sup.cst.sub.DN] < 1 does not necessarily imply that [[rho].sup.var.sub.DN] < 1. In other words, the fact that the algorithm with constant coefficients (the interface values) theoretically converges does not ensure that the algorithm with variable coefficients will. Indeed,

(2.12)[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

2.3. Convergence factor of the Robin-Robin algorithm with spatially variable coefficients. In this section we determine the convergence factor [[rho].sup.var.sub.RR] for the more general case of Robin-Robin interface conditions.

THEOREM 2.2 (Convergence factor with Robin-Robin transmission conditions). The convergence factor [[rho].sup.var.sub.RR] of the Schwarz algorithm (1.2) with Robin-Robin transmission conditions is

(2.13) [[rho].sup.var.sub.RR] = [absolute value of [([[lambda].sub.1] + [[lambda].sub.2])[K.sub.1] - 1] [([[lambda].sub.1] + [[lambda].sub.2])[K.sub.2] - 1]],

with [[lambda].sub.j] being the Fourier symbol of the operator [[LAMBDA].sub.j] in (1.2) and

(2.14) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and the eigenfunctions [[PHI].sub.n,j] and eigenvalues [c.sub.n,j] are solutions of the Sturm-Liouville problem (2.4).

Proof. Thanks to (2.6) and (2.7), we can express [[??].sub.1] and [[??].sub.2] in a compact form for the iterate k as

(2.15) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [K.sub.j] is given in (2.14). The problem on the interface x = 0 is given by the relations

(2.16) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and by combining (2.15) and (2.16), we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

By substituting those expressions in (2.16), we finally get a relation linking [??] and [??]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which leads to an expression for the convergence factor

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We note that this expression of the convergence factor is consistent with the expression (1.3) obtained in the case of constant (but discontinuous) coefficients. Indeed, if we set [[??].sub.1](x) = [[??].sub.2] (x) =0 in (2.13), we then have [K.sub.j] = 1/[square root of (i[omega][D.sub.j](0), which leads to (1.3) because [D.sub.j][[sigma].sup.[+ or -].sub.j] = [+ or -][square root of i[omega][D.sub.j](0)]. A convenient formulation of [[rho].sup.var.sub.RR], without complex numbers can be found in Appendix A. To conclude this section, we look at the asymptotic beha ior of [[rho].sup.var.sub.RR] and we can easily find that

[[rho].sup.var.sub.RR] ([omega] [right arrow] 0) = [[rho].sup.var.sub.RR] ([omega] [right arrow] [infinity]) [right arrow] 1,

which shows that the effect of the variability of the diffusion coefficients asymptotically vanishes when a Robin-Robin algorithm is used.

3. Numerical results. In this section we verify numerically the validity of the theoretical results presented in Section 2. To do this, we first briefly describe the rationale for the spatial variability of the diffusion coefficient and provide a typical profile which we will use for the numerical tests. Then we design a few experiments to illustrate the relevance of our theoretical results.

3.1. Planetary boundary layer turbulence. Unlike boundary layers in many engineering flows, the atmospheric and oceanic planetary boundary layers are almost always turbulent and cannot be explicitly resolved due to the insufficient vertical resolution in computational models. The numerical representation of those layers thus relies on the Reynolds decomposition: the flow is split into a mean (resolved) part <u> and a fluctuating (subgrid) part u' (where u can either represent a velocity component or an active tracer). When this decomposition is applied to nonlinear (advective) terms, this gives rise to additional terms and hence to a closure problem. The dominant expression in the turbulent boundary layers arising from the Reynolds decomposition is the divergence of the vertical part of (u 'w ') (where w denotes the vertical component of the velocity). Typically, this turbulent vertical flux is expressed as a function of the mean (resolved) part of the flow by using the down-gradient assumption, (u 'w ') = -D(x)[[partial derivative].sub.x] (u), where D(x) is the so-called eddy diffusivity or eddy-viscosity if u represents a velocity. This assumption explains why a one-dimensional diffusion equation, like the one studied in the present paper, is generally sufficient to locally represent the turbulent mixing in the boundary layers. The eddy diffusivity D(x) is defined to allow the flow to make the transition between its surface (the air-sea interface) and its interior (below the boundary layer) properties. This is the reason why D(x) exhibits a strong spatial variability. In this context, several ways to specify the coefficient D(x) have been proposed. The formulation most commonly used in the state-of-the-art numerical models can be found in [7] and [14]. Those formulations define the eddy diffusivity as

(3.1)[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with [h.sub.bl] the thickness of the boundary layer (depending on the state of the flow) and A a parameter setting the intensity of the mixing (note that D(x) is continuous and differentiable at x = [h.sub.bl]). Throughout this section we assume that D(x) is given by (3.1), and a typical profile for A = 0.5 m/s and [h.sub.bl] = 150 m is given in Figure 3.1.

[FIGURE 3.1 OMITTED]

In the remainder of this section we study first a Dirichlet-Neumann algorithm and then a Robin-Robin algorithm. We define spatially variable coefficients with AT =0.1 m/s (respectively [A.sub.2] = 0.5 m/s) and [h.sub.bl,1] = 50 m (respectively [h.sub.bl,2] = 150 m) on [[OMEGA].sub.1] (respectively [[OMEGA].sub.2]). The values of [v.sub.1] and [v.sub.2] (corresponding to the surface values [D.sub.1](0) and [D.sub.2] (0)) are chosen to be the same as the values used in [9] in the constant coefficient case. If we introduce [gamma] = [square root of ([[v.sub.2]/[v.sub.1]])], we investigate the two cases [gamma] =10 and [gamma] =[square root of ([square root of (10)]] with [v.sub.2] = 0.5 [m.sup.2]/s (the value of [v.sup.1] is adjusted depending on the value of [gamma]). Those various parameter values lead to diffusion profiles that can be found in the atmospheric and oceanic boundary layers. The discretization of the problem, the computational grid, as well as the initial conditions are described in [9, Section 5]. We use [DELTA]t = 100 s and a random initial guess on the interface so that it contains a wide range of the temporal frequencies that can be resolved by the computational grid.

3.2. Test case 1: Dirichlet-Neumann. The analytical convergence factor [[rho].sup.var.sub.DN] ([omega]) in Equation (2.8) is shown for different values of [gamma] in Figure 3.2. The eigenvalues [c.sub.n] and eigenfunctions [[PHI].sub.n] are computed numerically on the same computational grid as the model problem. We remark that depending on the jump in the coefficients through the interface, the spatial variability of the diffusivities either tend to accelerate the convergence speed (for [gamma] = [square root of ([square root of (10)]]) or to slow it down (for [gamma] = 10) compared to the convergence speed obtained with constant coefficients. As expected, the convergence factor for spatially variable coefficients is no longer independent of [omega], and for low frequencies we get a significant departure from the convergence rate of the algorithm with constant coefficients. The trend seen in the convergence factor [[rho].sup.var.sub.DN] ([omega]) determined at a continuous level is confirmed by the numerical results (Figure 3.2, top). These results as well as the asymptotic expression (2.12) call for caution when we use a Dirichlet-Neumann algorithm with spatially variable coefficients because it can lead to significantly different performances compared to the one obtained with constant coefficients. We can expect a Robin-Robin type algorithm to provide a more robust alternative thanks to the tuning of the [[lambda].sub.j] parameters.

[FIGURE 3.2 OMITTED]

3.3. Test case 2: Robin-Robin. Here, we denote by [[lambda].sup.cst.sub.j] the optimal Robin parameters obtained using the analytical results found in [9] for constant coefficients. We consider that these constant coefficients are the interface values [D.sub.j](0). Moreover, we denote by [[lambda].sup.var.sub.j] the Robin parameters optimized by solving numerically the problem (1.4) with the convergence factor [[rho].sup.var.sub.RR] as given in (2.13). This optimization is done using the Rosenbrock method [13] and by taking the parameters [[lambda].sup.cst.sub.j] to initialize the algorithm. We see from Figure 3.3 that the use of the parameters [[lambda].sup.var.sub.j] provide slightly better convergence properties compared to the parameters [[lambda].sup.cst.sub.j], regardless of the value of [gamma]. As for the Dirichlet-Neumann algorithm, we can check that our analytical study at the continuous level provides a convergence factor [[rho].sup.var.sub.RR] representative of the behavior of the algorithm at a discrete level (Figure 3.3, bottom). From Figure 3.3 (bottom left), we also see that the way we initialize the algorithm (with a random initial guess for [u.sup.0.sub.2] (0, t), t [member of] [0,T]) leads to the generation of a large range of temporal frequencies and more particularly low frequencies slowing down the convergence speed of the simulation using the parameters [[lambda].sup.cst.sub.j], although the latter provide a faster convergence than the parameters [[lambda].sup.var.sub.j] for most of the frequency spectrum. For our model problem, the use of the parameters [[lambda].sup.var.sub.j] provides a relatively modest improvement over the [[lambda].sup.cst.sub.j] parameters. However, in general, this statement has to be mitigated because if we consider [h.sub.bl,2] = 10 m instead of [h.sub.bl,2] = 150 m, we see in Figure 3.4 that the parameters obtained through an optimization of [[rho].sup.var.sub.RR] are clearly superior to the parameters [[lambda].sup.cst.sub.j]. For this case, we show in Figure 3.5 the asymptotic behavior of the optimized convergence rate and the associated Robin parameters [[lambda].sup.var.sub.j]. We numerically show that thanks to the theoretical work presented earlier in the paper, we are able to bring the proper adjustments to the Robin parameters so that our algorithm asymptotically maintains a good efficiency even in the presence of spatially variable coefficients.

[FIGURE 3.3 OMITTED]

4. Conclusion. In this paper we present and analyze a new approach to study the convergence properties of a global-in-time Schwarz algorithm in the case of a one-dimensional diffusion problem with spatially variable diffusion coefficients. We analytically derive an expression for the evolution of the errors of such an algorithm with respect to the iterates. Thanks to our formulation, we are able to gain a better understanding of the behavior of the associated convergence factor. We exhibit some interesting features that were not shown by the usual convergence studies with constant diffusion coefficients. We put particular emphasis on the fact that for low temporal frequencies, it can be a very inaccurate assumption to replace a variable diffusion coefficient by its constant interface value in the convergence study. Moreover, we also show that depending on the type of algorithm under consideration (Dirichlet-Neumann or Robin-Robin) the variability of the coefficients may have more or less impact on the asymptotic convergence properties. To be more attractive for practical applications, our approach requires further developments by performing an accurate study of the eigenvalue problems to improve our knowledge of the behavior of these eigenvalues with respect to the perturbations of the diffusion profiles.

[FIGURE 3.4 OMITTED]

[FIGURE 3.5 OMITTED]

Appendix A. Determination of the convergence factor in the case of variable coefficients. We recall (2.13):

(A.1) [rho] = [absolute value of [([[lambda].sub.1] + [lambda].sub.2])[K.sub.1] - 1] [([lambda].sub.1] + [lambda].sub.2])[K.sub.2] - 1]],

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

(A.1) can be rewritten as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

In order to determine the real and imaginary parts of [K.sub.1] and [K.sub.2], we can decompose each term appearing in the preceding expressions:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thanks to these equalities, we can recast Kj into the following form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and by noting that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we obtain

[K.sub.1] = ([e.sub.1](1 - [g.sub.1]) + [f.sub.1][h.sub.1]) + i([f.sub.1](1 - [g.sub.1]) - [e.sub.1][h.sub.1]),

[K.sub.2] = ([e.sub.2](1 - [g.sub.2]) + [f.sub.1][h.sub.2]) + i([f.sub.2](1 - [g.sub.2]) - [e.sub.1][h.sub.2]).

Hence, thanks to (A.2), this yields an expression for the convergence factor p without complex numbers.

Acknowledgments. This research was partially supported by the ANR project "COMMA" and by the INRIA project-team MOISE. The authors are thankful for the comments and suggestions of two anonymous reviewers, which helped to improve the clarity of this manuscript.

REFERENCES

[1] P. B. BAILEY, W. N. EVERITT, AND A. ZETTL, Regular and singular Sturm-Liouville problems with coupled boundary conditions, Proc. Roy. Soc. Edinburgh Sect. A, 126 (1996), pp. 505-514.

[2] V. W. EKMAN, On the influence of'the earth's rotation on ocean currents, Ark. Mat. Astron. Fys., 2 (1905), pp. 1-53.

[3] B. ENGQUIST AND A. MAJDA, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp., 31 (1977), pp. 629-651.

[4] M.J. GANDER ANDL. HALPERN, Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems, SIAM J. Numer. Anal., 45 (2007), pp. 666-697.

[5] M. J. GANDER, L. HALPERN, AND F. NATAF, Optimal convergence for overlapping and non-overlapping Schwarz waveform relaxation, in Eleventh International Conference on Domain Decomposition Methods London 1998, C. H. Lai, P. E. Bjorstad, M. Cross, and O. Widlund, eds., DDM.org, Augsburg, 1999, pp. 27-36.

[6] Q. KONG AND A. ZETTL, Eigenvalues ofregular Sturm-Liouville problems, J. Differential Equations, 131 (1996), pp. 1-19.

[7] W. G. LARGE, J. C. MCWILLIAMS, AND S. C. DONEY, Oceanic vertical mixing: A review and a model with a nonlocal boundarylayer parameterization, Rev. Geophys., 32 (1994), pp. 363-403.

[8] F. LEMARIE, Algorithmes de Schwarz et couplage ocean-atmosphere, Ph.D. Thesis, Laboratoire Jean Kuntz mann, Joseph Fourier University, Grenoble, 2008.

[9] F. LEMARlE, L. DEBREU, AND E. BLAYO, Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 1: the constant coefficients case, Electron. Trans. Numer. Anal., 40 (2013), pp. 148-169. http://etna.math.kent.edu/vol.40.2 013/pp14 8-16 9.dir

[10] P.-L. LIONS, On the Schwarz alternating method. III. A variant for nonoverlapping subdomains, in Third International Symposium on Domain Decomposition Methods for Partial Differential Equations, T. Chan, R. Glowinski, J. Periaux, and O. Widlund, eds., SIAM, Philadelphia, 1990, pp. 202-223.

[11] F. NATAF, F. ROGIER, AND E. DE STURLER, Optimal interface conditions for domain decomposition methods, Internal Report 301, CMAP, Ecole Polytechnique, Palaiseau, 1994.

[12] A. QUARTERONI AND A. VALLI, Domain Decomposition Methods for Partial Differential Equations, Numerical Mathematics and Scientific Computation, Oxford University Press, New York, 1999.

[13] H. H. ROSENBROCK, An automatic method for finding the greatest or least value of a function, Comput. J., 3 (1960), pp. 175-184.

[14] I. B. TROEN AND L. MAHRT, A simple model of the atmospheric boundary layer; sensitivity to surface evaporation, Boundary-Layer Meteorol., 37 (1986), pp. 129-148.

FLORIAN LEMARIE[dagger], LAURENT DEBREU[dagger], AND ERIC BLAYO[double dagger]

* Received September 29, 2010. Accepted March 11, 2013. Published online on July 5, 2013. Recommended by Martin Gander.

[dagger] INRIA Grenoble Rhone-Alpes, Montbonnot, 38334 Saint Ismier Cedex, France and Jean Kuntzmann Laboratory, BP 53, 38041 Grenoble Cedex 9, France, ({florian. lemarie, laurent. debreu}@inria .fr).

[double dagger] University of Grenoble, Jean Kuntzmann Laboratory, BP 53, 38041 Grenoble Cedex 9, France (eric.blayo@imag.fr).

Printer friendly Cite/link Email Feedback | |

Author: | Lemarie, Florian; Debreu, Laurent; Blayo, Eric |
---|---|

Publication: | Electronic Transactions on Numerical Analysis |

Article Type: | Report |

Date: | Jan 1, 2013 |

Words: | 5982 |

Previous Article: | Toward an optimized global-in-time Schwarz algorithm for diffusion equations with discontinuous and spatially variable coefficients. Part 1: the... |

Next Article: | Verified stability analysis using the Lyapunov matrix equation. |

Topics: |