# Optimized Robin Transmission Conditions for Anisotropic Diffusion on Arbitrary Meshes

Stella Krell, Martin J. Gander, Laurence Halpern, Florence Hubert

► To cite this version:

Stella Krell, Martin J. Gander, Laurence Halpern, Florence Hubert. Optimized Robin Transmission Conditions for Anisotropic Diffusion on Arbitrary Meshes. Domain Decomposition Methods in Science and Engineering XXVII, Jul 2022, Prague, Czech Republic. hal-05310505

HAL Id: hal-05310505

https://hal.science/hal-05310505v1

Submitted on 12 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimized Robin Transmission Conditions for Anisotropic Diffusion on Arbitrary Meshes

Martin J. Gander, Laurence Halpern, Florence Hubert, and Stella Krell

## 1 Introduction

We are interested in solving in parallel anisotropic diffusion problems of the form

$$ \mathcal{L}u := -\operatorname{div}(A\nabla u) + \eta u = g \quad \text{in } \Omega \subset \mathbb{R}^2, \quad u = 0 \quad \text{on } \partial\Omega, \quad (1) $$

where A is a symmetric positive definite matrix with $W^{1,\infty}$ coefficients,

$$ (x, y) \in \Omega \mapsto A(x, y) = \begin{pmatrix} A_{xx} & A_{xy} \\ A_{xy} & A_{yy} \end{pmatrix}, $$

and $(x, y) \in \Omega \mapsto \eta(x, y) \ge 0$ is in $L^\infty(\Omega)$. Schwarz algorithms for such problems are naturally formulated and studied at the continuous level. For a decomposition of the domain $\Omega$ into possibly non-overlapping subdomains $\Omega_j$, $j = 1, 2, \dots, J$, the parallel optimized Schwarz algorithm with Robin transmission conditions for the anisotropic diffusion problem (1) computes for iteration index $\ell = 1, 2, \dots$

$$ \begin{array}{ll} \mathcal{L}u_j^\ell = g & \text{in } \Omega_j, \\ u_j^\ell = 0 & \text{on } \partial\Omega_j \cap \partial\Omega, \\ A\nabla u_j^\ell \cdot \mathbf{n}_j + p u_j^\ell = -A\nabla u_i^{\ell-1} \cdot \mathbf{n}_i + p u_i^{\ell-1} & \text{on } \Gamma_{ji}, \end{array} \quad (2) $$

Martin J. Gander

Section de Mathématiques, Université de Genève, e-mail: martin.gander@unige.ch

Laurence Halpern

LAGA, Université Sorbonne Paris-Nord, e-mail: halpern@math.univ-paris13.fr

Florence Hubert

Aix-Marseille Université, CNRS, e-mail: florence.hubert@univ-amu.fr

Stella Krell

Université Côte d'Azur, Inria, CNRS, LJAD, e-mail: stella.krell@univ-cotedazur.fr
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Three typical discretizations for two subdomain decompositions: square-square (ss), triangle-square (ts) and triangle-quadrangle (tq).</figcaption></figure>

where $\mathbf{n}_j$ denotes the unit outer normal on the boundary of $\Omega_j$, and $\Gamma_{ji}$ denotes the portion of the interface where $\Omega_j$ takes data from $\Omega_i$. The efficiency of the algorithm is known to depend on the choice of the parameter $p$, which is usually optimized for a simple two subdomain decomposition, see [3] for the Laplace case. In [5], we showed at the continuous level for a general constant diffusion matrix $A$ that for $\Omega := (-a, a) \times (0, b)$ decomposed into two non-overlapping subdomains $\Omega_1 := (-a, 0) \times (0, b)$ and $\Omega_2 := (0, a) \times (0, b)$ with the interface $\Gamma_{12} = \Gamma_{21} := \partial\Omega_1 \cap \partial\Omega_2$, the optimized parameters and associated convergence factors are of the form

$$p^* = \sqrt{\tilde{f}(k_{\min})\tilde{f}(k_{\max})}, \quad \rho^* = \frac{\sqrt{\tilde{f}(k_{\max})} - \sqrt{\tilde{f}(k_{\min})}}{\sqrt{\tilde{f}(k_{\max})} + \sqrt{\tilde{f}(k_{\min})}}, \quad (3)$$

where for a general constant diffusion matrix $A$

$$\tilde{f}(k) := f(r(k)) \quad \text{with} \quad r(k) := \frac{1}{A_{xx}} \sqrt{\eta A_{xx} + \left(\frac{\pi k}{b}\right)^2 \det A}, \quad (4)$$

with the function $f$ defined for unbounded and bounded domains by

$$f(r) := \begin{cases} f_{\infty}(r) := A_{xx}r & a = \infty, \\ f_a(r) := f_{\infty}(r) \coth(ar) & a < \infty. \end{cases} \quad (5)$$

For both cases, the smallest frequency is $k_{\min} = 1$ and the largest frequency can be estimated by $k_{\max} = \frac{b}{h_y}$ for cell centered (cc) discretization, and $k_{\max} = \frac{b}{h_y} - 1$ for vertex centered (vc) discretizations, which are almost the same for small mesh size $h_y$ in the y direction, see below for more information.

We show for the three example meshes in Figure 1 the numerically computed convergence factors $\tilde{\rho}$ in Table 1 when running the optimized Schwarz algorithm discretized by Discrete Duality Finite Volumes (DDFV, see [5] for the DDFV Schwarz algorithm, and [7, 2, 1] for DDFV discretizations in general) for the Laplace problem, $A(x, y) = I$, and four anisotropic diffusion matrices, and characteristic mesh size $h_x = h_y =: h = \frac{1}{16}$, i.e the meshes in Figure 1 twice refined. We used the theoretically optimized value $p^* = p^*_{\infty, cvc}$ from (3) with $k_{\max} = \frac{b}{h_y} - 1$ corresponding to the vc scheme (index cvc for continuous vertex centered), see the comment at the end of section 3, and then also determined the numerically best working parameter $\tilde{\rho}^*$ and associated convergence factor $\tilde{\rho}^*$, which we computed (throughout the paper) per-
<table><caption>**Table 1** Numerically measured convergence factors $\tilde{\rho}$ of the optimized Schwarz algorithm for the three example meshes square-square (ss), triangle-square (ts) and triangle-quadrangle (tq) for the Laplace problem and four anisotropic diffusion problems with the theoretical parameter $p_{\infty,cvc}^*$ and the numerically best working one $\tilde{\rho}^*$.</caption><thead><tr><th colspan="2">Problem</th><th colspan="3">ss</th><th colspan="3">ts</th><th colspan="3">tq</th></tr><tr><th>A<sub>xx</sub></th><th>A<sub>yy</sub></th><th>p<sup>*</sup><sub>&infin;,cvc</sub></th><th>p&#773;</th><th>p&#773;</th><th>p&#773;</th><th>p&#773;</th><th>p&#773;</th><th>p&#773;</th><th>p&#773;</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>12.87</td><td>0.592</td><td>0.592</td><td>0.593</td><td>11.89</td><td>0.567</td><td>10.87</td><td>0.566</td><td>11.63</td><td>0.559</td></tr><tr><td>16</td><td>1</td><td>51.50</td><td>0.452</td><td>0.521</td><td>0.602</td><td>49.84</td><td>0.439</td><td>46.29</td><td>0.475</td><td>44.79</td><td>0.556</td></tr><tr><td>16</td><td>1/16</td><td>16.01</td><td>0.351</td><td>0.343</td><td>0.586</td><td>23.50</td><td>0.174</td><td>19.88</td><td>0.254</td><td>11.07</td><td>0.487</td></tr><tr><td>1</td><td>16</td><td>50.35</td><td>0.821</td><td>0.744</td><td>0.687</td><td>75.14</td><td>0.732</td><td>57.22</td><td>0.712</td><td>57.61</td><td>0.647</td></tr><tr><td>1/16</td><td>16</td><td>12.59</td><td>0.949</td><td>0.919</td><td>0.891</td><td>26.84</td><td>0.884</td><td>22.46</td><td>0.841</td><td>21.52</td><td>0.842</td></tr></tbody></table>

forming each time 100 iterations and using the last 40 to fit the linear convergence, to avoid initial fluctuations due to starting with a random initial guess.

We see from this experiment that for the Laplace problem the theoretically determined best parameter at the continuous level $p_{\infty,cvc}^*$ performs very well on all meshes, and is close to the numerically best working one $\tilde{\rho}^*$, with $\tilde{\rho} \approx \tilde{\rho}^*$. For anisotropic diffusion however this is not the case: the performance now depends on the mesh structure, and the numerically optimized parameter $\tilde{\rho}^*$ can be rather different from the theoretical parameter $p_{\infty,cvc}^*$. It is this difference we want to better understand, in particular for DDFV discretizations, which are highly accurate for anisotropic diffusion.

To start with our investigation, we plot in Figure 2 an example subdomain solution on the right subdomain $\Omega_2$ with interface value equal to 1 and vanishing source term for the Laplace case and two anisotropic diffusion cases. We see that the anisotropy deforms the solution quite a bit, and for $A_{xx}$ large, the subdomain clearly sees the boundary conditions at the outer boundary $\partial\Omega$ (Figure 2 middle), whereas for $A_{yy}$ large a boundary layer is forming close to the interface $\Gamma_{21}$ (Figure 2 right). This indicates that both the subdomain size, as well as the discretization, i.e. the mesh size, should influence the behavior of the optimized Schwarz method for anisotropic diffusion, and thus the best value of the parameter $p$.

<figure><img src="image_3.png" /><figcaption>**Fig. 2** Solutions for $A_{xx} = 1$, $A_{yy} = 1$ (left), for $A_{xx} = 16$, $A_{yy} = 1$ (middle) for $A_{xx} = 1$, $A_{yy} = 16$ (right), on an isotropic mesh.</figcaption></figure>
## 2 Optimized parameters at the discrete level

For rectangular meshes and for a diagonal anisotropy ($A_{xy} = 0$), it is easy to see (see e.g. [4]) that the DDFV scheme leads to two decoupled classical finite difference schemes, a cell centered (cc) scheme with unknowns at the cell centers, and a vertex centered (vc) scheme with unknowns at the vertices. In [4], we performed the optimization analysis in the same rectangular domain configuration as above, for a discretization associated to the step sizes $h_x$ and $h_y$ for both the cc and vc schemes for unbounded ($a = \infty$) and bounded ($a < \infty$) domains. The optimized parameters and associated convergence factors are again of form (3), with

$$ \tilde{f}(k) := f(v(k)), \quad v(k) := -\ln(\lambda(k)), \quad \lambda(k) := 1 + \frac{\mu(k)}{2} - \sqrt{\mu(k) + \frac{\mu(k)^2}{4}}, \quad (6) \\ \mu(k) := \frac{h_x^2}{A_{xx}} \left( 4 \frac{A_{yy}}{h_y^2} \sin^2 \left( \frac{k \pi h_y}{2b} \right) + \eta \right), $$

and the function $f$ is defined for the cc and vc schemes on unbounded and bounded domains by

$$ f(v) := \begin{cases} f_{\infty,cc}(v) := 2 \frac{A_{xx}}{h_x} \tanh \left( \frac{v}{2} \right), & a = \infty, \\ f_{a,cc}(v) := f_{\infty,cc}(v) \coth \left( \frac{av}{h_x} \right), & a < \infty, \\ f_{\infty,vc}(v) := \frac{A_{xx}}{h_x} \sinh (v), & a = \infty, \\ f_{a,vc}(v) = f_{\infty,vc}(v) \coth \left( \frac{av}{h_x} \right), & a < \infty. \end{cases} \quad (7) $$

Again the smallest frequency $k_{\min} = 1$, and the maximum frequencies can be estimated by $k_{\max} = \frac{b}{h_y}$ for the cc scheme and $k_{\max} = \frac{b}{h_y} - 1$ for the vc scheme.

## 3 Asymptotic analysis

In order to understand the difference in the performance of the optimized Schwarz method in the anisotropic case, we now present a new asymptotic analysis of the optimized parameters and associated convergence factors. We look at the asymptotic behavior as $h_x$ and $h_y$ tend to zero, their ratio being constant.

We start with the asymptotic analysis of the optimization results (4)–(5) at the continuous level. When inserting the smallest frequency $k = k_{\min}$ into (4)–(5), we get in the unbounded domain case

$$ \tilde{f}_{\infty}(k_{\min}) = \sqrt{\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A}, $$

and in the bounded domain case
$$ \tilde{f}_a(k_{\min}) = \sqrt{\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A} \coth \left( \frac{a}{A_{xx}} \sqrt{\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A} \right). $$

At the largest frequency $k = k_{\max}$, we obtain the same asymptotics, namely

$$ \tilde{f}_{\infty}(k_{\max}) = \tilde{f}_a(k_{\max}) = \frac{\pi \sqrt{\det A}}{h_y} + O(1). \quad (8) $$

Now when $h_y$ tends to zero, we see from (4) that $ar(k_{\max})$ tends to infinity, and therefore $coth(ar(k_{\max})) = 1 + o(h_y)$. We thus obtain for the unbounded domain case $a = \infty$ for the optimized parameter and associated convergence factor

$$ \begin{aligned} p_{\infty}^* &\sim \left(\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A\right)^{\frac{1}{4}} \left(\pi \sqrt{\det A}\right)^{\frac{1}{2}} h_y^{-\frac{1}{2}}, \\ \rho_{\infty}^* &\sim 1 - 2 \left(\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A\right)^{\frac{1}{4}} \left(\pi \sqrt{\det A}\right)^{-\frac{1}{2}} h_y^{\frac{1}{2}}, \end{aligned} $$

where $f(h_y) \sim g(h_y)$ means $\lim_{h_y \to 0} \frac{f(h_y)}{g(h_y)} = 1$, and when $a < \infty$, we get

$$ \begin{aligned} p_a^* &\sim \left(\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A\right)^{\frac{1}{4}} \left(\pi \sqrt{\det A}\right)^{\frac{1}{2}} \left(\coth\left(\frac{a}{A_{xx}} \sqrt{\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A}\right)\right)^{\frac{1}{2}} h_y^{-\frac{1}{2}}, \\ \rho_a^* &\sim 1 - 2 \left(\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A\right)^{\frac{1}{4}} \left(\pi \sqrt{\det A}\right)^{-\frac{1}{2}} \left(\coth\left(\frac{a}{A_{xx}} \sqrt{\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 \det A}\right)\right)^{\frac{1}{2}} h_y^{\frac{1}{2}}. \end{aligned} $$

We see that the asymptotic behavior in the mesh size is the same, but the constants differ between the bounded and unbounded domain case, clearly indicating that the continuous analysis on the bounded domain can take into account the anisotropy observed in Figure 2.

We next perform an asymptotic analysis of the optimization results (6) and (7) at the discrete level. For a diagonal diffusion matrix $A$, at the minimum frequency, $k = k_{\min}$, we obtain from (6)

$$ \mu(k_{\min}) = \frac{h_x^2}{A_{xx}} \left( 4 \frac{A_{yy}}{h_y^2} \sin^2\left(\frac{\pi h_y}{2b}\right) + \eta \right) = \frac{h_x^2}{A_{xx}^2} \left( \eta A_{xx} + \left(\frac{\pi}{b}\right)^2 A_{xx} A_{yy} + O(h_y^2) \right). $$

Hence $\mu(k_{\min}) \to 0$ when the mesh is refined, and because $\lambda(k_{\min}) \sim 1 - \sqrt{\mu(k_{\min})}$ and $\tilde{f}(k_{\min}) \sim \frac{A_{xx}}{h_x} \sqrt{\mu(k_{\min})}$, we obtain

$$ \tilde{f}_{\infty,cc}(k_{\min}) \sim \tilde{f}_{\infty,vc}(k_{\min}) \sim \sqrt{\eta A_{xx} + \left(\frac{\pi}{b}\right)^2 A_{xx} A_{yy}}. \quad (9) $$

At the highest frequency, $k = k_{\max}$, we obtain for the cc scheme
$$ \mu_{cc}(k_{max}) = \frac{h_x^2}{A_{xx}} \left( 4 \frac{A_{yy}}{h_y^2} \sin^2\left(\frac{\pi}{2}\right) + \eta \right) = \frac{h_x^2}{A_{xx}^2} \left( \eta A_{xx} + 4 \frac{A_{xx} A_{yy}}{h_y^2} \right) \sim 4\beta, $$

where $\beta := \frac{A_{yy} h_x^2}{h_y^2 A_{xx}}$, and similarly for the vc scheme,

$$ \mu_{vc}(k_{max}) = \frac{h_x^2}{A_{xx}} \left( 4 \frac{A_{yy}}{h_y^2} \sin^2\left(\frac{\pi}{2}(1-h_y)\right) + \eta \right) = \frac{h_x^2}{A_{xx}^2} \left( 4 \frac{A_{xx} A_{yy}}{h_y^2} + O(1) \right) \sim 4\beta. $$

Note that the case of a Laplacian with an isotropic square mesh corresponds to the parameter value $\beta = 1$. By hyperbolic trigonometric calculus, and $\frac{A_{xx}}{h_x} = \frac{\sqrt{A_{xx}A_{yy}}}{h_y\sqrt{\beta}}$, we obtain the alternative formula $f_{\infty,cc}(\nu(k)) = 2\frac{A_{xx}}{h_x} \frac{1-\lambda(k)}{1+\lambda(k)}$, which yields

$$ \begin{aligned} \tilde{f}_{\infty,cc}(k_{max}) &= 2 \frac{A_{xx}}{h_x} \frac{-\beta + \sqrt{\beta + \beta^2}}{1 + \beta - \sqrt{\beta + \beta^2}} = \frac{\sqrt{A_{xx}A_{yy}}}{h_y \sqrt{\beta}} 2 \frac{-\beta + \sqrt{\beta + \beta^2}}{1 + \beta - \sqrt{\beta + \beta^2}} \\ &= \frac{\sqrt{A_{xx}A_{yy}}}{h_y \sqrt{\beta}} 2 \frac{\sqrt{\beta + \beta^2}}{1 + \beta} := \frac{\sqrt{A_{xx}A_{yy}}}{h_y} \psi_{cc}(\beta), \end{aligned} $$

with $\psi_{cc}(\beta) = \frac{2}{\sqrt{1+\beta}}$. Similarly, since $f_{\infty,vc}(\nu(k)) = \frac{A_{xx}}{h_x} \frac{1-\lambda(k)^2}{2\lambda(k)}$ by hyperbolic trigonometric calculus, we obtain

$$ \begin{aligned} f_{\infty,vc}(k_{max}) &= \frac{A_{xx}}{h_x} \frac{2(-\beta + \sqrt{\beta + \beta^2})(1 + \beta - \sqrt{\beta + \beta^2})}{1 + 2\beta - 2\sqrt{\beta + \beta^2}} \\ &= \frac{\sqrt{A_{xx}A_{yy}}}{h_y \sqrt{\beta}} 2(-\beta + \sqrt{\beta + \beta^2})(1 + \beta - \sqrt{\beta + \beta^2})(1 + 2\beta + 2\sqrt{\beta + \beta^2}) \\ &= \frac{\sqrt{A_{xx}A_{yy}}}{h_y \sqrt{\beta}} 2\sqrt{\beta + \beta^2} := \frac{\sqrt{A_{xx}A_{yy}}}{h_y} \psi_{vc}(\beta), \end{aligned} $$

with $\psi_{vc}(\beta) = 2\sqrt{1+\beta}$. Note that in the special case $\beta = 1$, we get $\psi_{cc}(\beta) = \sqrt{2}$ and $\psi_{vc}(\beta) = 2\sqrt{2}$, a factor 2 difference. For the unbounded domain case, $a = \infty$, we then obtain for the optimized parameters and associated convergence factors of the cc and vc schemes

$$ p_{\infty,cc}^* \sim \psi_{cc}(\beta)^{\frac{1}{2}} \sqrt{A_{xx}A_{yy}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} h_y^{-\frac{1}{2}}, $$

$$ p_{\infty,vc}^* \sim \psi_{vc}(\beta)^{\frac{1}{2}} \sqrt{A_{xx}A_{yy}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} h_y^{-\frac{1}{2}}, $$

$$ \rho_{\infty,cc}^* \sim 1 - 2\psi_{cc}(\beta)^{-\frac{1}{2}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} h_y^{\frac{1}{2}}, $$

$$ \rho_{\infty,vc}^* \sim 1 - 2\psi_{vc}(\beta)^{-\frac{1}{2}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} h_y^{\frac{1}{2}}. $$
In the bounded domain case, $a < \infty$, we see that $\coth(\frac{av(k_{\max})}{h_x}) \sim 1$ and when $\mu(k_{\min}) \to 0$, we have $v(k_{\min}) \sim -\sqrt{\mu(k_{\min})}$, which implies

$$ \frac{av(k_{\min})}{h_x} \sim \frac{a}{\sqrt{A_{xx}}} \sqrt{\eta + \left(\frac{\pi}{b}\right)^2 A_{yy}} \Rightarrow \coth\left(\frac{av(k_{\min})}{h_x}\right) \sim \coth\left(\frac{a}{\sqrt{A_{xx}}} \sqrt{\eta + \left(\frac{\pi}{b}\right)^2 A_{yy}}\right). \quad (10) $$

We therefore get for the optimized parameters and associated convergence factors for the cc and vc schemes in the bounded domain case

$$ p_{a,cc}^* \sim \psi_{cc}(\beta)^{\frac{1}{2}} \sqrt{A_{xx}A_{yy}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} \coth \left( \frac{a}{\sqrt{A_{xx}}} \sqrt{\eta + \left(\frac{\pi}{b}\right)^2 A_{yy}} \right)^{\frac{1}{2}} h_y^{-\frac{1}{2}}, $$

$$ p_{a,vc}^* \sim \psi_{vc}(\beta)^{\frac{1}{2}} \sqrt{A_{xx}A_{yy}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} \coth \left( \frac{a}{\sqrt{A_{xx}}} \sqrt{\eta + \left(\frac{\pi}{b}\right)^2 A_{yy}} \right)^{\frac{1}{2}} h_y^{-\frac{1}{2}}, $$

$$ \rho_{a,cc}^* \sim 1 - 2\psi_{cc}(\beta)^{-\frac{1}{2}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} \coth \left( \frac{a}{\sqrt{A_{xx}}} \sqrt{\eta + \left(\frac{\pi}{b}\right)^2 A_{yy}} \right)^{\frac{1}{2}} h_y^{\frac{1}{2}}, $$

$$ \rho_{a,vc}^* \sim 1 - 2\psi_{vc}(\beta)^{-\frac{1}{2}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} \coth \left( \frac{a}{\sqrt{A_{xx}}} \sqrt{\eta + \left(\frac{\pi}{b}\right)^2 A_{yy}} \right)^{\frac{1}{2}} h_y^{\frac{1}{2}}. $$

These formulas take both the domain size and the mesh resolution into account, also when the mesh is not chosen appropriately for the anisotropy under consideration.

If one can not use separate parameters for the cc and vc components in a DDFV implementation, it was shown in [4] that the optimized choice for one parameter is of the form

$$ p_{a,ddfv}^* = \sqrt{f_{a,cc}(v(k_{\min}))f_{a,vc}(v(k_{\max}))}, $$

and since asymptotically we have $f_{a,cc}(v(k_{\min})) \sim f_{a,vc}(v(k_{\min}))$ from (9) and (10), one should use the optimized parameter $p_{a,vc}^* \sim p_{a,ddfv}^*$ in that case.

The continuous and discrete asymptotic results lead to the following general theorem.

### Theorem 1 (Optimized Robin parameter for diagonal anisotropic diffusion)

The optimized Schwarz method (2) for the anisotropic diffusion problem (1) with diagonal diffusion matrix A and a subdomain decomposition of the rectangle $\Omega = (-a, a) \times (0, b)$ into two non-overlapping subdomains $\Omega_1 := (-a, 0) \times (0, b)$ and $\Omega_2 := (0, a) \times (0, b)$ has for small mesh size $h_y$ the asymptotically optimized parameter and associated convergence factor

$$ p^* \sim \psi^{\frac{1}{2}} \sqrt{A_{xx}A_{yy}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} c^{\frac{1}{2}} h_y^{-\frac{1}{2}}, \quad (11) $$

$$ \rho^* \sim 1 - 2\psi^{-\frac{1}{2}} \left( \frac{\eta}{A_{yy}} + \left(\frac{\pi}{b}\right)^2 \right)^{\frac{1}{4}} c^{\frac{1}{2}} h_y^{\frac{1}{2}}, \quad (12) $$

where in the unbounded domain case, $a = \infty$, we have $c = 1$, whereas in the bounded domain case, $a < \infty$, we have
<figure><img src="image_4.png" /><figcaption>Fig. 3 Graph of the functions $\psi_{cc}(\beta)$ and $\psi_{vc}(\beta)$ for the discrete analysis, compared to $\psi(\beta) = \pi$ (dotted) from the continuous analysis for small and large $\beta$ range.</figcaption></figure>

$$c := c(a, b, A_{xx}, A_{yy}, \eta) = \coth \left( \frac{a}{\sqrt{A_{xx}}} \sqrt{\eta + \left(\frac{\pi}{b}\right)^2 A_{yy}} \right). \quad (13)$$

Furthermore, in the continuous case $\psi = \pi$, and in the discrete case we have

$$\psi := \psi_{cc}(\beta) = \frac{2}{\sqrt{1 + \beta}} \quad \text{or} \quad \psi = \psi_{vc}(\beta) := 2\sqrt{1 + \beta} \quad (14)$$

for the cell centered or vertex centered discretizations, with

$$\beta := \frac{A_{yy} h_x^2}{h_y^2 A_{xx}}. \quad (15)$$

Plotting the $\psi(\beta)$ functions in Figure 3, we see that if $\beta = 1$ then the continuous and discrete analyses give about the same optimized parameter $p^*$ and associated convergence factor, especially for the vc scheme. Since $\beta = \frac{A_{yy} h_x^2}{h_y^2 A_{xx}}$, this can be achieved by having equal mesh sizes $h_x = h_y$ and isotropic diffusion $A_{xx} = A_{yy}$, or by adapting the mesh sizes to the anisotropy, $h_y^2 = \frac{A_{xx}}{A_{yy}} h_x^2$. Such an adaptation is also recommended for accuracy, since a Taylor expansion gives

$$A_{xx} \frac{u(x+h_x,y)-2u(x,y)-u(x-h_x,y)}{h_x^2} + A_{yy} \frac{u(x,y+h_y)-2u(x,y)-u(x,y-h_y)}{h_y^2} = (A_{xx}\partial_{xx} + A_{yy}\partial_{yy})u(x,y) + \frac{1}{12}(A_{xx}h_x^2\partial_x^4 + A_{yy}h_y^2\partial_y^4)u(x,y) + \dots, \quad (16)$$

and from the separation of variables solution $u(x, y) = e^{-\frac{k\pi}{b}\sqrt{\frac{A_{yy}}{A_{xx}}x}} \sin(\frac{k\pi}{b}y)$ we see that the fourth derivative in x scales like $\frac{A_{yy}^2}{A_{xx}^2}$, while the fourth derivative in y does not scale in these entries, and hence to balance the error term, we should choose
$$A_{xx}h_x^2 \frac{A_{yy}^2}{A_{xx}^2} \approx A_{yy}h_y^2 \quad \implies \quad \frac{h_x^2}{A_{xx}} \frac{A_{yy}}{h_y^2} = \beta \approx 1. \qquad (17)$$

Hence for $\beta \approx 1$, we can use the continuous analysis results and expect good performance, also in highly anisotropic cases, provided the mesh is adapted accordingly. If $\beta$ is very different from one, we should use the parameters from the discrete analysis to get good performance. We also see from Figure 3 (right) that for large $\beta$ the optimized parameters for the cc and vc schemes are becoming more and more different, and (12) together with (14) indicates that the cc scheme is converging much faster than the vc scheme in these not well resolved mesh situations. In the DDFV case with general meshes, where both cc and vc discretizations are involved, the importance will then lie on a good optimization of the vc parameter, the cc parameter playing only a secondary role in these not well resolved cases.

Next, we see from Theorem 1 that if $c \approx 1$, then we can use the unbounded domain analysis, since the only term depending on the domain bound $a$ on the left and right is $c$. Now $c \approx 1$ if the argument of the coth is large, i.e. either the domains and thus $a/b$ is large, or $\eta$ is large, or $A_{yy}/A_{xx}$ is large, which is illustrated in Figure 1 on the right, where we see that the outer boundary on the right does not play a major role any more<fnref n="1" />. If none of these hold, then the bounded domain analysis needs to be used to obtain good performance.

Finally, from $\rho^*$ in Theorem 1, we see the algorithm will converge very fast with the well chosen $p^*$, provided $A_{yy}$ is small or $\eta$ large, or $\psi(\beta)$ is small. Having $\psi(\beta)$ small is however not advisable, because the discretization accuracy is only good for $\beta \approx 1$, see (17).

# 4 Numerical experiments

We can now explain the discrepancies we observed in Table 1 as soon as we solve anisotropic diffusion problems. There are two reasons: the first one is that when using the optimized parameter $p^*$ from the continuous, unbounded domain analysis, the fact that the subdomains are actually bounded in a concrete computation becomes important as soon as the diffusion in the orthogonal direction to the interface is large, and the cross diffusion tangential to the interface is small. This is visible also in Figure 1 showing a corresponding solution in the middle, where we can clearly see that the boundary on the right makes the solution decay linearly in the direction orthogonal to the interface, in stark contrast to the Laplace case on the left in Figure 1, where the decay is exponential. The second reason for discrepancies is the uniform discretization, which can not resolve well the boundary layers close to the top and bottom boundaries in Figure 1 (middle), and close to the interface in Figure 1 (right) which also influence the convergence of the Schwarz method.

<footnote n="1">¹ For example, in the case $A_{xx} = 1$ and $A_{yy} = 16$, the difference is of order $10^{-11}$.</footnote>
<table><caption>**Table 2** Results corresponding to Table 1 but now using the theoretical parameter $p_{a,cvc}^*$ from the bounded domain analysis.</caption><thead><tr><th colspan="2">Problem</th><th colspan="3">ss</th><th colspan="3">ts</th><th colspan="3">tq</th></tr><tr><th>A<sub>xx</sub></th><th>A<sub>yy</sub></th><th>P<sup>*</sup><sub>a,cvc</sub></th><th>ρ̃</th><th>ρ̃</th><th>ρ̃</th><th>ρ̃<sup>*</sup></th><th>ρ̃<sup>*</sup></th><th>ρ̃<sup>*</sup></th><th>ρ̃<sup>*</sup></th><th>ρ̃<sup>*</sup></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>12.48</td><td>0.582</td><td>0.581</td><td>0.583</td><td>11.89</td><td>0.567</td><td>10.87</td><td>0.566</td><td>11.63</td><td>0.559</td></tr><tr><td>16</td><td>1</td><td>60.59</td><td>0.514</td><td>0.578</td><td>0.651</td><td>49.84</td><td>0.439</td><td>46.29</td><td>0.475</td><td>44.79</td><td>0.556</td></tr><tr><td>16</td><td>1/16</td><td>28.04</td><td>0.258</td><td>0.436</td><td>0.741</td><td>23.50</td><td>0.174</td><td>19.88</td><td>0.254</td><td>11.07</td><td>0.487</td></tr><tr><td>1</td><td>16</td><td>48.75</td><td>0.826</td><td>0.751</td><td>0.695</td><td>75.14</td><td>0.732</td><td>57.22</td><td>0.712</td><td>57.61</td><td>0.647</td></tr><tr><td>1/16</td><td>16</td><td>12.19</td><td>0.950</td><td>0.921</td><td>0.894</td><td>26.84</td><td>0.884</td><td>22.46</td><td>0.841</td><td>21.52</td><td>0.842</td></tr></tbody></table>

<table><caption>**Table 3** Results corresponding to Table 2 but now using the discrete theoretical parameters $p_{a,cc}^*$ and $p_{a,vc}^*$, and the numerically best working ones $\tilde{p}_{cc}^*$ and $\tilde{p}_{vc}^*$.</caption><thead><tr><th colspan="2">Problem</th><th colspan="3">ss</th><th colspan="3">ts</th><th colspan="3">tq</th></tr><tr><th>A<sub>xx</sub></th><th>A<sub>yy</sub></th><th>P<sup>*</sup><sub>a,cc</sub></th><th>P<sup>*</sup><sub>a,vc</sub></th><th>ρ̃</th><th>ρ̃</th><th>ρ̃</th><th>ρ̃<sup>*</sup></th><th>ρ̃<sup>*</sup></th><th>ρ̃<sup>*</sup></th><th>ρ̃<sup>*</sup></th></tr></thead><tbody><tr><td>1</td><td>1</td><td>8.62</td><td>12.22</td><td>0.573</td><td>0.572</td><td>0.574</td><td>8.62</td><td>11.93</td><td>0.566</td><td>7.73</td><td>11.38</td><td>0.533</td><td>10.49</td><td>10.49</td><td>0.527</td></tr><tr><td>16</td><td>1</td><td>49.16</td><td>50.56</td><td>0.444</td><td>0.509</td><td>0.592</td><td>49.59</td><td>49.87</td><td>0.439</td><td>45.87</td><td>45.89</td><td>0.468</td><td>39.61</td><td>40.13</td><td>0.514</td></tr><tr><td>16</td><td>1/16</td><td>23.48</td><td>23.48</td><td>0.174</td><td>0.347</td><td>0.698</td><td>23.50</td><td>23.44</td><td>0.173</td><td>19.75</td><td>20.24</td><td>0.242</td><td>11.42</td><td>11.65</td><td>0.466</td></tr><tr><td>1</td><td>16</td><td>19.07</td><td>84.09</td><td>0.723</td><td>0.728</td><td>0.733</td><td>20.01</td><td>80.71</td><td>0.714</td><td>44.46</td><td>66.21</td><td>0.653</td><td>13.78</td><td>58.50</td><td>0.621</td></tr><tr><td>1/16</td><td>16</td><td>1.84</td><td>54.59</td><td>0.806</td><td>0.834</td><td>0.861</td><td>1.13</td><td>51.09</td><td>0.796</td><td>1.90</td><td>36.72</td><td>0.756</td><td>0.69</td><td>30.80</td><td>0.733</td></tr></tbody></table>

As a first remedy, we use the optimized parameter $p_a^*$ from the continuous, bounded domain analysis to take into account the boundedness of the domains. From Table 2 we see that this already improves the performance of the method when the diffusion is large in the orthogonal direction to the interface and small tangentially. However for the other cases using the bounded domain analysis is not sufficient due to the bad mesh resolution in the anisotropic case.

We therefore now use the discrete optimized formulas $p_{a,vc}^*$ and $p_{a,cc}^*$ in our DDFV Schwarz code, which are perfectly adapted to the anisotropy of the problem we are solving on bounded subdomains, and truly optimize both the vc and cc scheme component convergence also for the not well chosen mesh resolution. We show the corresponding results in Table 3. We see that now our parameters predicted by the discrete analysis for the cc and vc schemes give performance close to the truly best possible ones for rectangular meshes, and still work well on general meshes for which our analysis is not valid any more. Furthermore, the performance still follows our asymptotic analysis, as the plots of the convergence factors under mesh refinement in Figure 4 indicate.

We finally show numerical results using an anisotropic mesh which gives better approximate discrete solutions, see the truncation error analysis in (16). We show the corresponding results for such meshes in Table 4, and in Figure 5. We see that the continuous analysis gives now very good predictions for the optimized parameters for the vc scheme, while for the cc scheme their value is still a bit overestimated. This does however not influence the performance very much.
<figure><img src="image_5.png" /><figcaption>Fig. 4 Asymptotic dependence of $1 - \tilde{\rho}$ on the mesh size for isotropic meshes and the anisotropic diffusion problems in Table 3, with $h = h_y = h_x$. From top left to bottom right: $(A_{xx}, A_{yy}) = (16, 1)$, $(16, \frac{1}{16})$, $(1, 16)$, $(\frac{1}{16}, 16)$.</figcaption></figure>

## 5 Conclusions

Using asymptotic analysis, we explained rigorously numerical observations on the performance of DDFV optimized Schwarz methods applied to anisotropic diffusion. We showed that for strong anisotropic diffusion solved on uniform, non-adapted meshes, one needs optimized parameters from a more subtle discrete analysis, continuous optimization does not suffice. When using suitably adapted, anisotropic meshes such that the discrete solution is a good approximation of the continuous one, optimized parameters from a continuous analysis perform however well. We also showed numerically that this remains true if one uses meshes for which a detailed asymptotic analysis as ours on Cartesian meshes can not be performed. For extensions of the DDFV Schwarz algorithm to Navier-Stokes problems, see [6].

<table><caption>Table 4 Results obtained using the discrete optimized parameters for adapted anisotropic meshes.</caption><thead><tr><th colspan="2">Problem</th><th colspan="8">ss aniso</th></tr><tr><th>A<sub>xx</sub></th><th>A<sub>yy</sub></th><th>P<sup>*</sup><sub>a,ccc</sub></th><th>P<sup>*</sup><sub>a,cvc</sub></th><th>P&#773;<sub>c</sub></th><th>P<sup>*</sup><sub>a,cc</sub></th><th>P<sup>*</sup><sub>a,vc</sub></th><th>P&#773;</th><th>P&#773;<sup>*</sup><sub>cc</sub></th><th>P&#773;<sup>*</sup><sub>vc</sub></th><th>P&#773;<sup>*</sup></th></tr></thead><tbody><tr><td>16</td><td>1</td><td>125.13</td><td>124.15</td><td>0.730</td><td>83.94</td><td>118.73</td><td>0.718</td><td>82.30</td><td>111.96</td><td>0.705</td></tr><tr><td>16</td><td><sup>1</sup>&frasl;<sub>16</sub></td><td>115.32</td><td>115.09</td><td>0.749</td><td>77.37</td><td>109.43</td><td>0.737</td><td>77.37</td><td>102.45</td><td>0.724</td></tr><tr><td>1</td><td>16</td><td>50.35</td><td>48.75</td><td>0.601</td><td>33.67</td><td>47.67</td><td>0.581</td><td>33.37</td><td>46.43</td><td>0.573</td></tr><tr><td><sup>1</sup>&frasl;<sub>16</sub></td><td>16</td><td>12.59</td><td>12.19</td><td>0.601</td><td>8.42</td><td>11.92</td><td>0.580</td><td>8.42</td><td>11.63</td><td>0.574</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Asymptotic dependence of $1 - \tilde{\rho}$ on the mesh size for anisotropic meshes and the anisotropic diffusion problems in Table 4. From top left to bottom right: ($A_{xx}, A_{yy}$) = (16, 1), (16, $\frac{1}{16}$), (1, 16), ($\frac{1}{16}$, 16).</figcaption></figure>

References

1. Andreianov, B., Boyer, F., and Hubert, F. Discrete duality finite volume schemes for Leray-Lions type problems on general 2D meshes. *Numerical Methods for PDE 23*(1), 145–195 (2007).

2. Domelevo, K. and Omnes, P. A finite volume method for the Laplace equation on almost arbitrary two-dimensional grids. *M2AN Math. Model. Numer. Anal.* **39**(6), 1203–1249 (2005).

3. Gander, M. J. Optimized Schwarz methods. *SIAM Journal on Numerical Analysis* **44**(2), 699–731 (2006).

4. Gander, M. J., Halpern, L., Hubert, F., and Krell, S. Discrete optimization of Robin transmission conditions for anisotropic diffusion with discrete duality finite volume methods. *Vietnam Journal of Mathematics* **49**, 1349–1378 (2021).

5. Gander, M. J., Halpern, L., Hubert, F., and Krell, S. Optimized Schwarz Methods for Anisotropic Diffusion with Discrete Duality Finite Volume Discretizations. *Moroccan Journal of Pure and Applied Analysis* **7**(2), 182–213 (2021).

6. Goudon, T., Krell, S., and Lissoni, G. Non-overlapping Schwarz algorithms for the incompressible Navier–Stokes equations with DDFV discretizations. *ESAIM: Mathematical Modelling and Numerical Analysis* **55**(4), 1271–1321 (2021).

7. Hermeline, F. Approximation of diffusion operators with discontinuous tensor coefficients on distorted meshes. *Comput. Methods Appl. Mech. Engrg.* **192**(16-18), 1939–1959 (2003).