Decoupling PDE Computation with Intrinsic or Inertial Robin Interface Condition

TL;DR

This paper introduces a novel intrinsic or inertial Robin interface condition for decoupling multi-domain, multi-physics PDE problems, providing a mathematically justified framework that improves numerical decoupling efficiency.

Contribution

It derives a new interface condition for PDE decoupling that is mathematically and physically justified, enhancing the effectiveness of numerical methods for multi-model problems.

Findings

The new Robin condition improves decoupling accuracy.

Numerical experiments confirm the effectiveness of the approach.

The framework generalizes decoupling methods for multi-physics applications.

Abstract

We study decoupled numerical methods for multi-domain, multi-physics applications. By investigating various stages of numerical approximation and decoupling and tracking how the information is transmitted across the interface for a typical multi-modeling model problem, we derive an approximate intrinsic or inertial type Robin condition for its semi-discrete model. This new interface condition is justified both mathematically and physically in contrast to the classical Robin interface condition conventionally introduced for decoupling multi-modeling problems. Based on the intrinsic or inertial Robin condition, an equivalent semi-discrete model is introduced, which provides a general framework for devising effective decoupled numerical methods. Numerical experiments also confirm the effectiveness of this new decoupling approach.

Figures

Tables

Table 1. Table 1: Errors of ‖ u h , N − u e ​ x ​ t ​ ( T ) ‖ 0 , Ω subscript norm subscript 𝑢 ℎ 𝑁 subscript 𝑢 𝑒 𝑥 𝑡 𝑇 0 Ω \|u_{h,N}-u_{ext}(T)\|_{0,\Omega} with ρ 1 = ρ 2 = 1 , β 1 = β 2 = β ​ ( x , y ) , Δ ​ t = h 2 . formulae-sequence subscript 𝜌 1 subscript 𝜌 2 1 subscript 𝛽 1 subscript 𝛽 2 𝛽 𝑥 𝑦 Δ 𝑡 superscript ℎ 2 \rho_{1}=\rho_{2}=1,~{}\beta_{1}=\beta_{2}=\beta(x,y),~{}\Delta t=h^{2}.

h	Coupled scheme	Decoupled iRN scheme
$\frac{1}{8}$	3.43107e-2	3.75366e-2
$\frac{1}{16}$	9.20988e-3	1.13574e-2
$\frac{1}{32}$	2.34374e-3	2.61407e-3
$\frac{1}{64}$	5.88545e-4	5.77793e-4

Table 2. Table 2: Errors of ‖ u h , N − u e ​ x ​ t ​ ( T ) ‖ 0 , Ω subscript norm subscript 𝑢 ℎ 𝑁 subscript 𝑢 𝑒 𝑥 𝑡 𝑇 0 Ω \|u_{h,N}-u_{ext}(T)\|_{0,\Omega} with ρ 1 = ρ 2 = 1 , β 1 = β 2 = β ​ ( x , y ) , Δ ​ t = h 2 . formulae-sequence subscript 𝜌 1 subscript 𝜌 2 1 subscript 𝛽 1 subscript 𝛽 2 𝛽 𝑥 𝑦 Δ 𝑡 superscript ℎ 2 \rho_{1}=\rho_{2}=1,~{}\beta_{1}=\beta_{2}=\beta(x,y),~{}\Delta t=h^{2}.

h	Coupled scheme	Decoupled DN scheme
$\frac{1}{8}$	3.43107e-2	3.53891e-2
$\frac{1}{16}$	9.20988e-3	9.46749e-3
$\frac{1}{32}$	2.34374e-3	2.40475e-3
h	Decoupled iRR scheme	Decoupled RR scheme
		${\alpha }_1=10,{\alpha }_2=5$	${\alpha }_1=1,{\alpha }_2=1$
$\frac{1}{8}$	3.53436e-2	3.24908e-2	3.80028e-2
$\frac{1}{16}$	9.55398e-3	8.81404e-3	1.39637e-2
$\frac{1}{32}$	2.25541e-3	2.26911e-3	5.08893e-3

Table 3. Table 3: Errors of ‖ u h , N − u e ​ x ​ t ​ ( T ) ‖ 0 , Ω subscript norm subscript 𝑢 ℎ 𝑁 subscript 𝑢 𝑒 𝑥 𝑡 𝑇 0 Ω \|u_{h,N}-u_{ext}(T)\|_{0,\Omega} with ρ 1 = 10 , ρ 2 = 1 , β 1 = β 2 = β ​ ( x , y ) , Δ ​ t = h 2 formulae-sequence formulae-sequence subscript 𝜌 1 10 formulae-sequence subscript 𝜌 2 1 subscript 𝛽 1 subscript 𝛽 2 𝛽 𝑥 𝑦 Δ 𝑡 superscript ℎ 2 \rho_{1}=10,~{}\rho_{2}=1,~{}\beta_{1}=\beta_{2}=\beta(x,y),~{}\Delta t=h^{2} .

h	Coupled scheme	Decoupled DN scheme
$\frac{1}{8}$	3.36012e-2	$\infty$
$\frac{1}{16}$	9.01136e-3	$\infty$
$\frac{1}{32}$	2.29268e-3	$\infty$
h	Decoupled iRR scheme	Decoupled RR scheme
		${\alpha }_1=10,{\alpha }_2=5$	${\alpha }_1=1,{\alpha }_2=1$
$\frac{1}{8}$	3.24630e-2	3.18627e-2	3.72954e-2
$\frac{1}{16}$	8.80706e-3	8.64307e-3	1.37006e-2
$\frac{1}{32}$	2.27580e-3	2.26911e-3	5.00114e-3

Table 4. Table 4: Errors of ‖ u h , N − u e ​ x ​ t ​ ( T ) ‖ 0 , Ω subscript norm subscript 𝑢 ℎ 𝑁 subscript 𝑢 𝑒 𝑥 𝑡 𝑇 0 Ω \|u_{h,N}-u_{ext}(T)\|_{0,\Omega} with ρ 1 = 1 , ρ 2 = 10 , β 1 = β 2 = β ​ ( x , y ) , Δ ​ t = h 2 . formulae-sequence formulae-sequence subscript 𝜌 1 1 formulae-sequence subscript 𝜌 2 10 subscript 𝛽 1 subscript 𝛽 2 𝛽 𝑥 𝑦 Δ 𝑡 superscript ℎ 2 \rho_{1}=1,~{}\rho_{2}=10,~{}\beta_{1}=\beta_{2}=\beta(x,y),~{}\Delta t=h^{2}.

h	Coupled scheme	Decoupled DN scheme
$\frac{1}{8}$	3.39101e-2	3.21780e-2
$\frac{1}{16}$	9.09826e-3	9.38034e-3
$\frac{1}{32}$	2.31506e-3	2.38207e-3
h	Decoupled iRR scheme	Decoupled RR scheme
		${\alpha }_1=10,{\alpha }_2=5$	${\alpha }_1=1,{\alpha }_2=1$
$\frac{1}{8}$	3.23888e-2	3.21780e-2	3.76394e-2
$\frac{1}{16}$	8.82536e-3	8.72945e-3	1.38221e-2
$\frac{1}{32}$	2.28889e-3	2.24812e-3	5.04232e-3

Equations

$${\rho }_1\frac{\partial u_1}{\partial t}-\nabla \cdot ({\beta }_1\nabla u_1)=f_1,$$

$${\rho }_2\frac{\partial u_2}{\partial t}-\nabla \cdot ({\beta }_2\nabla u_2)=f_2,$$

$$u_1=u_2,$$

$${\beta }_1\nabla u_1\cdot {\mathbf{n}}_1+{\beta }_2\nabla u_2\cdot {\mathbf{n}}_2=q,$$

$$u_1=u_{1D},$$

$$u_2=u_{2D},$$

$${\mathbf{V}}_1=\left\{v_1\in {\mathbf{H}}^1({\Omega }_1):v_1=0\text{mbox}on\partial {\Omega }_1\backslash \Gamma \right\},$$

$${\mathbf{V}}_2=\left\{v_2\in {\mathbf{H}}^1({\Omega }_2):v_2=0\text{mbox}on\partial {\Omega }_2\backslash \Gamma \right\},$$

$$\mathbf{V}=\left\{(v_1,v_2)\in {\mathbf{V}}_1\times {\mathbf{V}}_2:v_1=v_2\text{mbox}on\Gamma \right\},$$

$$\tilde{\mathbf{V}}=\left\{(v_1,v_2)\in {\mathbf{V}}_1\times {\mathbf{V}}_2\right\},$$

$${\mathbf{V}}^{\Gamma }=\left\{(v_1,v_2)\in {\mathbf{V}}_1\times {\mathbf{V}}_2:v_1=v_2=0\text{mbox}on\Gamma \right\},$$

$${\mathbf{\Lambda }}_{\Gamma }=\{\xi =v_1{|}_{\Gamma }:v_1\in {\mathbf{V}}_1\}=\{\xi =v_2{|}_{\Gamma }:v_2\in {\mathbf{V}}_2\},$$

$${\mathbf{V}}_1^{\Gamma }=\{v_1\in {\mathbf{V}}_1:v_1=0\text{mbox}on\Gamma \},$$

$${\mathbf{V}}_2^{\Gamma }=\{v_2\in {\mathbf{V}}_2:v_2=0\text{mbox}on\Gamma \}.$$

$$\begin{array}{r}{\rho }_1{\left(\frac{\partial u_1^C}{\partial t},v_1\right)}_{{\Omega }_1}+{\left({\beta }_1\nabla u_1^C,\nabla v_1\right)}_{{\Omega }_1}+{\rho }_2{\left(\frac{\partial u_2^C}{\partial t},v_2\right)}_{{\Omega }_2}+{\left({\beta }_2\nabla u_2^C,\nabla v_2\right)}_{{\Omega }_2}\\ ={\left(f_1,v_1\right)}_{{\Omega }_1}+{\left(f_2,v_2\right)}_{{\Omega }_2},\forall v=(v_1,v_2)\in \mathbf{V}.\end{array}$$

$$\begin{array}{r}{\rho }_1{\left(\frac{\partial u_{1,h}^S}{\partial t},v_{1,h}\right)}_{{\Omega }_1}+{\left({\beta }_1\nabla u_{1,h}^S,\nabla v_{1,h}\right)}_{{\Omega }_1}+{\rho }_2{\left(\frac{\partial u_{2,h}^S}{\partial t},v_{2,h}\right)}_{{\Omega }_2}+\\ {\left({\beta }_2\nabla u_{2,h}^S,\nabla v_{2,h}\right)}_{{\Omega }_2}={\left(f_{1,h},v_{1,h}\right)}_{{\Omega }_1}+{\left(f_{2,h},v_{2,h}\right)}_{{\Omega }_2},\forall v_h=(v_{1,h},v_{2,h})\in {\mathbf{V}}_h,\end{array}$$

$$A=\left[\begin{array}{ccccccc}{a}_{11}& a_{12}& & & & & \\ a_{21}& a_{22}& a_{23}& & & & \\ & \ddots & \ddots & \ddots & & & \\ & & a_{ii-1}& a_{ii}& a_{ii+1}& & \\ & & & \ddots & \ddots & \ddots & \\ & & & & \ddots & \ddots & \ddots \\ & & & & & a_{nn-1}& a_{nn}\end{array}\right]\stackrel{\text{Lumped-mass}}{\to }$$

$$A_L=\left[\begin{array}{ccccccc}{a}_1& 0& & & & & \\ 0& a_2& 0& & & & \\ & \ddots & \ddots & \ddots & & & \\ & & 0& a_i& 0& & \\ & & & \ddots & \ddots & \ddots & \\ & & & & \ddots & \ddots & \ddots \\ & & & & & 0& a_n\end{array}\right],$$

$$a_i=\sum _{j=1}^n{a}_{ij}.$$

$${\left(u_h,v_h\right)}_{\Omega }={\int }_{\Omega }{u}_h{v}_h,\text{for any}{u}_h,v_h\in W_h.$$

$$Q_{\tau }(f)=\frac{1}{3}\text{area}(\tau )\sum _{j=1}^{3}f(P_{\tau ,j}),$$

$${\left(u_h,v_h\right)}_{\Omega ,L}=\sum _{\tau \in {\mathcal{T}}_h}{Q}_{\tau }\left(u_h{v}_h\right).$$

$${\mathcal{L}}_{1,h}:{\mathbf{\Lambda }}_{\Gamma ,h}⟶{\mathbf{V}}_{1,h},$$

$${\mathcal{L}}_{2,h}:{\mathbf{\Lambda }}_{\Gamma ,h}⟶{\mathbf{V}}_{2,h},$$

$${\left({\mathbf{B}}_{1,h}{\xi }_h,{\lambda }_h\right)}_{\Gamma }={\left({\mathcal{L}}_{1,h}{\xi }_h,{\mathcal{L}}_{1,h}{\lambda }_h\right)}_{{\Omega }_1,L},\forall {\xi }_h,{\lambda }_h\in {\mathbf{\Lambda }}_{\Gamma ,h}$$

$${\left({\mathbf{B}}_{2,h}{\xi }_h,{\lambda }_h\right)}_{\Gamma }={\left({\mathcal{L}}_{2,h}{\xi }_h,{\mathcal{L}}_{2,h}{\lambda }_h\right)}_{{\Omega }_2,L},\forall {\xi }_h,{\lambda }_h\in {\mathbf{\Lambda }}_{\Gamma ,h}.$$

$${\left(v_{1,h},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1,L}={\left({\mathbf{B}}_{1,h}{v}_{1,h},{\xi }_h\right)}_{\Gamma }.$$

$${\left(v_{1,h},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1,L}={\left({\tilde{v}}_{1,h},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1,L}+{\left({\mathcal{L}}_{1,h}{v}_{1,h},{\mathcal{L}}_1{\xi }_h\right)}_{{\Omega }_1,L}={\left({\mathbf{B}}_{1,h}{v}_{1,h},{\xi }_h\right)}_{\Gamma }.$$

$${\left(v_{2,h},{\mathcal{L}}_{2,h}{\xi }_h\right)}_{{\Omega }_2,L}={\left({\mathbf{B}}_{2,h}{v}_{2,h},{\xi }_h\right)}_{\Gamma }.$$

$$\begin{array}{r}{\rho }_1{\left(\frac{\partial u_{1,h}^L}{\partial t},v_{1,h}\right)}_{{\Omega }_1,L}+{\left({\beta }_1\nabla u_{1,h}^L,\nabla v_{1,h}\right)}_{{\Omega }_1}+{\rho }_2{\left(\frac{\partial u_{2,h}^L}{\partial t},v_{2,h}\right)}_{{\Omega }_2,L}+\\ {\left({\beta }_2\nabla u_{2,h}^L,\nabla v_{2,h}\right)}_{{\Omega }_2}={\left(f_{1,h},v_{1,h}\right)}_{{\Omega }_1}+{\left(f_{2,h},v_{2,h}\right)}_{{\Omega }_2},\forall v_h=(v_{1,h},v_{2,h})\in {\mathbf{V}}_h.\end{array}$$

$$\Vert u_{i,h}^L-u_i^C{\Vert }_{L^2({\Omega }_i)}=O(h^2),i=1,2,$$

$$\Vert u_{i,h}^L-u_i^C{\Vert }_{H^1({\Omega }_i)}=O(h),i=1,2.$$

$$\begin{array}{r}{\rho }_1{\left(\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},v_{1,h}\right)}_{{\Omega }_1,L}+{\left({\beta }_1\nabla u_{1,h}^{\tilde{L}},\nabla v_{1,h}\right)}_{{\Omega }_1}+{\rho }_2{\left(\frac{\partial u_{2,h}^{\tilde{L}}}{\partial t},v_{2,h}\right)}_{{\Omega }_2,L}+\\ {\left({\beta }_2\nabla u_{2,h}^{\tilde{L}},\nabla v_{2,h}\right)}_{{\Omega }_2}={\left(f_{1,h},v_{1,h}\right)}_{{\Omega }_1}+{\left(f_{2,h},v_{2,h}\right)}_{{\Omega }_2},\forall v_h=(v_{1,h},v_{2,h})\in {\mathbf{V}}_h.\end{array}$$

$$\begin{array}{r}{\rho }_1{\left({\mathbf{B}}_{1,h}\frac{\partial u_{2,h}^{\tilde{L}}}{\partial t},{\xi }_h\right)}_{\Gamma }+\\ \left[{\rho }_2{\left(\frac{\partial u_{2,h}^{\tilde{L}}}{\partial t},{\mathcal{L}}_{2,h}{\xi }_h\right)}_{{\Omega }_2,L}+{\left({\beta }_2\nabla u_{2,h}^{\tilde{L}},\nabla ({\mathcal{L}}_{2,h}{\xi }_h)\right)}_{{\Omega }_2}-{\left(f_{2,h},{\mathcal{L}}_{2,h}{\xi }_h\right)}_{{\Omega }_2}\right]=\\ {\rho }_1{\left({\mathbf{B}}_{1,h}\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},{\xi }_h\right)}_{\Gamma }-\\ \left[{\rho }_1{\left(\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1,L}+{\left({\beta }_1\nabla u_{1,h}^{\tilde{L}},\nabla ({\mathcal{L}}_{1,h}{\xi }_h)\right)}_{{\Omega }_1}-{\left(f_{1,h},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1}\right]+\\ {\rho }_1{\left({\mathbf{B}}_{1,h}\frac{\partial (u_{1,h}^{\tilde{L}}-u_{2,h}^{\tilde{L}})}{\partial t},{\xi }_h\right)}_{\Gamma },\forall {\xi }_h\in {\mathbf{\Lambda }}_{\Gamma ,h}.\end{array}$$

$$\begin{array}{rl}& \left[{\rho }_1{\left(\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1,L}+{\left({\beta }_1\nabla u_{1,h}^{\tilde{L}},\nabla ({\mathcal{L}}_{1,h}{\xi }_h)\right)}_{{\Omega }_1}-{\left(f_{1,h},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1}\right]+\\ & \left[{\rho }_2{\left(\frac{\partial u_{2,h}^{\tilde{L}}}{\partial t},{\mathcal{L}}_{2,h}{\xi }_h\right)}_{{\Omega }_2,L}+{\left({\beta }_2\nabla u_{2,h}^{\tilde{L}},\nabla ({\mathcal{L}}_{2,h}{\xi }_h)\right)}_{{\Omega }_2}-{\left(f_{2,h},{\mathcal{L}}_{2,h}{\xi }_h\right)}_{{\Omega }_2}\right]=0.\end{array}$$

$${\rho }_1{\left(\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1,L}={\rho }_1{\left({\mathbf{B}}_{1,h}\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},{\xi }_h\right)}_{\Gamma }.$$

$${\rho }_1{\left(\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1,L}={\rho }_1{\left({\mathbf{B}}_{1,h}\frac{\partial u_{2,h}^{\tilde{L}}}{\partial t},{\xi }_h\right)}_{\Gamma }-{\rho }_1{\left({\mathbf{B}}_{1,h}\frac{\partial e_D}{\partial t},{\xi }_h\right)}_{\Gamma }.$$

$$\begin{array}{rl}& \left[{\rho }_1{\left({\mathbf{B}}_{1,h}\frac{\partial u_{2,h}^{\tilde{L}}}{\partial t},{\xi }_h\right)}_{\Gamma }+{\left({\beta }_1\nabla u_{1,h}^{\tilde{L}},\nabla ({\mathcal{L}}_{1,h}{\xi }_h)\right)}_{{\Omega }_1}-{\left(f_{1,h},{\mathcal{L}}_{1,h}{\xi }_h\right)}_{{\Omega }_1}\right]+\\ & \left[{\rho }_2{\left(\frac{\partial u_{2,h}^{\tilde{L}}}{\partial t},{\mathcal{L}}_{2,h}{\xi }_h\right)}_{{\Omega }_2,L}+{\left({\beta }_2\nabla u_{2,h}^{\tilde{L}},\nabla ({\mathcal{L}}_{2,h}{\xi }_h)\right)}_{{\Omega }_2}-{\left(f_{2,h},{\mathcal{L}}_{2,h}{\xi }_h\right)}_{{\Omega }_2}\right]=\\ & {\rho }_1{\left({\mathbf{B}}_{2,h}\frac{\partial e_D}{\partial t},{\xi }_h\right)}_{\Gamma }.\end{array}$$

Full text

{AMS}

65N30, 65M55, 74F10

1 Introduction

Multi-domain, multi-modeling problems find important applications in science and engineering. They usually involve coupled processes or systems in more than one geometric domains or physical fields, and arise from classical parallel computing, domain decomposition, to contemporary multi-scale and multi-physics computation. Mathematically, they are modeled by different partial differential equations(PDEs) in different domains with certain interface coupling conditions, such as fluid-structure interactions (FSI) [9, 16, 23], coupling fluid flow with porous media flow [1, 7, 10], etc..

Numerical methods for multi-domain, multi-modeling problems are generally classified as monolithic and partitioned approaches. The monolithic approach usually leads to coupled schemes [3, 20, 21, 24]. Although unconditionally stable and convergent in general, those coupled schemes usually result in significant difficulties and inflexibility in the design and choice of mesh generation, PDE discretization, algebraic solvers, as well as software development.

On the other hand, certain partitioned approaches, often called loosely coupled, explicit coupling, or decoupled schemes, have also been investigated in the literature [2, 4, 5, 6, 11, 12, 13, 17, 18], where the sub-model in each domain is solved locally so that existing legacy solvers may be applied by easy software integration. There are many other important physical and numerical considerations that appeal to partitioned approaches for treating different sub-models in their own physical regions independently, particularly because of the very different physical properties and time scales. However, decoupled approaches are usually theoretically much more difficult. Decoupling could easily lead to numerical instability or non-convergence. For instance, many decoupled methods in FSI applications are so-called unconditionally instable due to the artificial added-mass effect [8, 14]. Even with classical non-overlapping domain decomposition, for a long time people have been still searching for and trying to understand how to decouple the computation effectively to ensure stability and convergence. The most typical and conventional approach is to match the given interface conditions, for instance the continuity of the Dirichlet and Neumann data, in one way or another by solving local sub-models with certain given boundary conditions and then exchanging the information across the interface in certain ad hoc or intuitive way, which often leads to decoupled methods of the Dirichlet-Dirichlet or Dirichlet-Neumann type. These methods are sometimes stable and convergent, sometime not, or not very efficient. A Robin type condition has been introduced by combining the original Dirichlet and Neumann interface conditions to obtain a set of mathematically equivalent interface conditions, which results in some other types of decoupled methods such as Dirichlet-Robin, Neumann-Robin, and Robin-Robin type of schemes. Sometimes, this may help to improve the effectiveness and efficiency. However, such an approach is artificial and lack of physical justification, which results in both theoretical and practical numerical difficulties for the Robin type decoupling approach. For instance, in applications like FSI the Dirichlet data correspond to velocity, while the Neumann data correspond to interfacial force, therefore, there is a mismatch between these two different types of physical quantities. In the classical Robin approach, this physical mismatch is somehow meredied by mathematically introducing a relaxation parameter, which is to be tuned for the purposes of numerical stability and convergence. The optimal relaxation parameter is then compensated for the mathematical effect of the physical mismatch. However, the optimal parameter analysis is often very difficult, even for certain simple model problems. The choice and tuning of the relaxation parameter with the existing numerical methods are practically very difficult and even impossible for most of the real applications, never mentioning the theoretical analysis and understanding of the numerical stability and convergence.

In an interesting FSI application of a thin-wall structure coupled with a fluid flow, a so-call intrinsic Robin condition was introduced with certain formal derivation, which leads to a stable decoupled scheme and overcomes the added-mass instability effect [12, 13]. Motivated by this, we believe that this could provide certain insight and help understand the interface mechanism, as well as provide guidance for developing decoupling methodology for general multi-domain and multi-physics applications.

We investigate a typical multi-modeling model problem. By cruising various stages of numerical approximation and decoupling, and track how the information is transmitted across the interface, we clearly demonstrate how the interface conditions evolve and affect the data and error transmission across the interface during spacial discretization, time discretization, as well as decoupling. An approximate intrinsic or inertial type Robin condition is rigorously derived for a semi-discrete model with finite element spacial discretization, which makes sense both mathematically and physically in contrast to the classical Robin condition. Based on this new interface condition, an equivalent semi-discrete model is introduced, which provides a general framework for devising effective decoupled numerical methods. Numerical experiments also confirm the effectiveness of this new decoupling approach.

The paper is organized as follows. The multi-domain model problem is described in Section 2. The derivation and analysis of the new intrinsic or inertial type Robin interface condition are carried out in Section 3. The decoupling approach based on the intrinsic or inertial type Robin interface condition is introduced in Section 4, followed by numerical experiments in Section 5. Concluding remarks are given in Section 6.

2 Model Problem

We consider a scalar coupled model consisting of two heat equations in two adjacent subdomains coupled through interface conditions. Let $\Omega_{1}=(0,1)\times(0,1)$ and $\Omega_{2}=(1,2)\times(0,1)$ be the two subdomains, with the interface $\Gamma=\partial\Omega_{1}\cap\partial\Omega_{2}$ as showed in Figure 1, where ${\bm{n}}_{1}$ and ${\bm{n}}_{2}$ are the unit outward normal vectors to $\partial\Omega_{1}$ and $\partial\Omega_{2}$, respectively.

The coupled PDE model is described as:

[TABLE]

where $u_{1}$ and $u_{2}$ are temperature, $\rho_{1}$ and $\rho_{2}$ are density, $\beta_{1}$ and $\beta_{2}$ are material diffusivity, $f_{1}$ and $f_{2}$ are external source, and the Dirichlet and Neumann conditions are imposed on the interface by (3) and (4). For illustration and without loss of generality, we will assume $q=0$.

It serves as a typical model for investigating decoupling techniques for various important applications, bridging classical domain decomposition and parallel computing with contemporary multi-physics computing. When $\rho_{1}=\rho_{2}$ and $\beta_{1}=\beta_{2}$, this coupled model is equivalent to a underline single model of a heat equation defined over the global domain, being partitioned to subdomain problems with non-overlapping domain decomposition so that various domain decomposition time marching schemes or iterative methods may be devised for parabolic or elliptic problems with parallel computing. When $\beta_{1}\neq\beta_{2}$, and/or $\rho_{1}\neq\rho_{2}$, it may be used as a model problem to study decoupling techniques for jumping coefficients and multi-scale applications. The scalar model problem may also be used to investigate and understand the numerical and physical properties of decoupling techniques for being further extended to other important applications in multi-modeling, multi-physics, multi-scale computation such as FSI.

Define the function spaces

[TABLE]

Integration by parts gives the weak form of the coupled model (1)-(6):

**Model C: The Weak Form of the Continuous Coupled Model

**

Find $u^{C}=(u_{1}^{C},u_{2}^{C})\in{\bm{V}}$, such that

[TABLE]

3 Semi-Discretization and Approximate Intrinsic or Inertial Type Robin Conditions

The weak formulation is mathematically equivalent to the strong form of the original coupled PDE model in the conventional sense, and its solution still satisfies the original Dirichlet and Neumann conditions (2.3) and (2.4) on the interface. There could be various equivalent interface conditions to be derived for various purposes. To understand what should be the appropriate interface conditions for a numerical scheme to satisfy for the sake of stability and convergence, in particular when we want to devise effective decoupling techniques later, let us start with the original continuous model and investigate the mechanism on how the interface conditions would change from the original model due to discrete approximations.

Note that the issue of numerical stability during time marching in numerical computation is caused by the discretization in both space and time. To understand how the spacial discretizatin affects the interface conditions, yet without the stability factor, let us concentrate on the finite element semi-discrete approximation:

**Model $S_{h}$: The Semi-Discrete Model with Finite Element Approximation in Space

**

Find $u_{h}^{S}=(u_{1,h}^{S},u_{2,h}^{S})\in{\bm{V}}_{h}$, such that

[TABLE]

where ${\bm{V}}_{h}$ is a standard finite element subspace of ${\bm{V}}$. Note now that Model $S_{h}$ is no longer equivalent to Model C, but just an approximation. Yet there are standard error estimates for the solutions of these two models under certain regularity assumptions.

However, it is essentially important to point out that the approximate solution $u_{h}^{S}$ of Model $S_{h}$ no longer satisfies the interface conditions of the original coupled model although the Dirichlet condition remains enforced in the solution space ${\bm{V}}_{h}$, to which little attention has been paid in the past, mainly because the fully implicit time discretization will then lead to an unconditionally stable scheme, and the convergence and error analysis are ready in the classical finite element analysis, so it appears unnecessary to further investigate what happens on the interface. However, if decoupling is necessary as a further approximation, it would then become crucially important to find out what is indeed the corresponding NEW interface coupling conditions implied by this approximate semi-discrete model due to the spacial discretization, which would then provide us insight for decoupling the numerical computation properly in order to keep maintaining numerical stability and convergence for decoupled methods.

Unfortunately, it is difficult to derive the local interface conditions for such a weak formulation defined in the integration form over the entire geometric domain. For finding out such an interface condition approximately, motivated by [13], we propose to apply the lumped-mass technique [22], which will lead to a further approximation to the semi-discrete model.

Let us now briefly review the lumped-mass technique. Let $A$ be a $n\times n$ matrix and $A_{L}$ be the corresponding lumped-mass matrix

[TABLE]

where

[TABLE]

Denote $W_{h}\subset H^{1}$ consists of continuous, piecewise linear functions on a quasiuniform family of triangulations $\mathcal{T}_{h}$ of any domain $\Omega$, where $\Omega=\Omega_{1}$ or $\Omega=\Omega_{2}$. The inner product in $W_{h}$ is defined by

[TABLE]

Let $\tau$ be a triangle of the triangulation $\mathcal{T}_{h}$ and $P_{\tau,j}$ be its vertices. Define the quadrature formula for any function $f$

[TABLE]

then define the lumped-mass approximation of inner product $\left(u_{h},v_{h}\right)_{\Omega}$ by

[TABLE]

Let $\{v_{i}\}_{i\leq n}\in W_{h}$ be the finite element basis, the following two properties hold [22]:

• •

$\left(v_{i},v_{j}\right)_{\Omega,L}=0$ for $i\neq j$,

• •

$\left(v_{i},v_{i}\right)_{\Omega,L}=\sum_{j=1}^{n}\left(v_{i},v_{j}\right)_{\Omega}$.

From the properties, if $A$ is generated by $\left(v_{i},v_{j}\right)_{\Omega}$, then $A_{L}$ is generated by $\left(v_{i},v_{j}\right)_{\Omega,L}$.

Denote ${\bm{V}}_{1,h},{\bm{V}}_{1,h},{\bm{\Lambda}}_{\Gamma,h},{\bm{V}}_{1,h}^{\Gamma},{\bm{V}}_{2,h}^{\Gamma}$ be the corresponding discrete spaces with continuous piecewise linear functions. Applying the lumped-mass technique to ${\bm{V}}_{1,h},{\bm{V}}_{1,h}$, we may build connections among the inner products in ${\bm{\Lambda}}_{\Gamma,h}$ and ${\bm{V}}_{1,h},~{}{\bm{V}}_{2,h}$. Define two discrete lifting operators

[TABLE]

such that the nodal values of $\mathcal{L}_{1,h}\xi_{h}$ and $\mathcal{L}_{2,h}\xi_{h}$ vanish out of $\Gamma$ for all $\xi_{h}\in{\bm{\Lambda}}_{\Gamma,h}$. The interface operator ${\bm{B}}_{1,h}:~{}{\bm{\Lambda}}_{\Gamma,h}\longrightarrow{\bm{\Lambda}}_{\Gamma,h}$ is defined by ${\bm{B}}_{1,h}=\mathcal{L}_{1,h}^{\ast}\mathcal{L}_{1,h}$, where $\mathcal{L}_{1,h}^{\ast}$ is the adjoint operator of $\mathcal{L}_{1,h}$ w.r.t the lumped-mass inner product in $V_{1,h}$. ${\bm{B}}_{2,h}$ is also defined similarly. Based on the discrete operators introduced above, we have

[TABLE]

The following Lemmas are useful in the analysis below. For simplicity, we denote $\mathcal{L}_{1}(v_{1}|_{\Gamma}),~{}\mathcal{L}_{2}(v_{2}|_{\Gamma})$ by $\mathcal{L}_{1}(v_{1}),~{}\mathcal{L}_{2}(v_{2})$.

Lemma 3.1**.**

For all $\left(v_{1,h},\xi\right)\in{\bm{V}}_{1,h}\times{\bm{\Lambda}}_{\Gamma,h}$, we have

[TABLE]

Proof 3.2**.**

For $v_{1,h}\in{\bm{V}}_{1,h}$, it can be decomposed into $v_{1,h}=\tilde{v}_{1,h}+\mathcal{L}_{1,h}v_{1,h}$ with $\tilde{v}_{1,h}\in{\bm{V}}_{1,h}^{\Gamma}$. From the definition of lumped-mass inner product (13), we can conclude that $\left(\tilde{v}_{1,h},\mathcal{L}_{1,h}\xi\right)_{\Omega_{1},L}=0$ for all $\xi_{h}\in{\bm{\Lambda}}_{\Gamma,h}$ since all the nodal values of $\tilde{v}_{1,h}\cdot\mathcal{L}_{1,h}\xi_{h}$ are [math]. Therefore, using (15), we have

[TABLE]

Similarly, we also have

Lemma 3.3**.**

For all $\left(v_{2,h},\xi_{h}\right)\in{\bm{V}}_{2,h}\times{\bm{\Lambda}}_{\Gamma,h}$, we have

[TABLE]

We may now define the lumped-mass semi-discrete model as a further approximation, denoted by Model $L_{h}$.

**Model $L_{h}$: The Lumped-Mass Semi-Discrete Model

**

Find $u_{h}^{L}=(u_{1,h}^{L},u_{2,h}^{L})\in{\bm{V}}_{h}$, such that

[TABLE]

Under certain regularity assumptions, there are also standard error estimates for the solutions of the two semi-discrete models, as well as the original continuous model [22]. In particular, if $u^{C}=\left(u_{1}^{C},u_{2}^{C}\right)$ has $H^{2}$ regularity, for each subdomain we have

[TABLE]

There are many important applications in numerical computation where the Dirichlet condition enforced in the solution space ${\bm{V}}_{h}$ should be relaxed at certain steps for decoupling purposes. Such cases arise from matching the Dirichlet data on the interface by the classical non-overlapping domain decomposition approach for constructing parallel time marching schemes for parabolic problems or for constructing parallel iterative methods for elliptic problems by simulating the time evolution to the equilibrium state, to today’s decoupling for multi-physics computation. For understanding the interface mechanism in order to construct effective decoupled numerical methods as well as to conduct the corresponding numerical analysis in the future, it will be fundamentally important to examine the following modified problem to the above Model $L_{h}$ by removing the Dirichlet interface condition from the solution space ${\bm{V}}_{h}$, where the modification of ${\bm{V}}_{h}$ is defined as ${\tilde{\bm{V}}}_{h}=\left\{(v_{1,h},v_{2,h})\in{\bm{V}}_{1,h}\times{\bm{V}}_{2,h}\right\}$.

**Problem $P_{h}$: A Perturbation to Model $L_{h}$

**

Consider a function $u_{h}^{\tilde{L}}=(u_{1,h}^{\tilde{L}},u_{2,h}^{\tilde{L}})\in{\tilde{\bm{V}}}_{h}$ satisfying

[TABLE]

The following theorem is revealing.

Theorem 3.4**.**

If $u_{h}^{\tilde{L}}\in{\tilde{\bm{V}}}_{h}$ satisfies (20) in Problem $P_{h}$, then $u_{h}^{\tilde{L}}$ satisfies the following condition:

[TABLE]

Proof 3.5**.**

Given $\xi_{h}\in{\bm{\Lambda}}_{\Gamma,h}$, taking $v_{1,h}=\mathcal{L}_{1,h}\xi_{h}$ and $v_{2,h}=\mathcal{L}_{2,h}\xi_{h}$ in (20), we have

[TABLE]

The lumped-mass property (16) in Lemma 3.1 allows us to condense the inner product in subdomain $\Omega_{1}$ to the interface $\Gamma$:

[TABLE]

In order to pass the Dirichlet information from $\Omega_{1}$ to $\Omega_{2}$, denote $u_{1,h}^{\tilde{L}}=u_{2,h}^{\tilde{L}}-e_{D}$, we have

[TABLE]

Applying (24) to (22) gives

[TABLE]

Reorganizing the terms by separating the information on each subdomain yields

[TABLE]

Note that the left hand side in (26) consists of two part: the expression in the bracket corresponds to the terms in (20) from Problem $P_{h}$ for subdomain $\Omega_{2}$, and the time derivative term being condensed to the interface from subdomain $\Omega_{1}$, but then exchanged the Dirichlet interface data with subdomain $\Omega_{2}$ yet leaving an error term in the right hand side, which will play the key role in the intrinsic Robin condition to be derived.

For symmetry and later on to compare with the semi-discrete model, let us add and subtract the same term $\rho_{1}\left(\frac{\partial u_{1,h}^{\tilde{L}}}{\partial t},\mathcal{L}_{1,h}\xi_{h}\right)_{\Omega_{1},L}$ to equation (26), which gives

[TABLE]

*Using (23) again, the added inner product term outside of the bracket in the right hand side of the above equation may also be condensed to the interface, which yields (21) and thus completes the proof. *

The fundamental significance of $Theorem$ 3.4 as well as its derivation is in that the error equation (21) reveals the interface mechanism of the semi-discrete model, in particular the impact of the error in the Dirichlet data on the interface conditions when the Dirichlet data is transmitted across the interface, which results instead in the error of the time derivative of the Dirichlet data in the form of $\rho_{1}{\bm{B}}_{1,h}\frac{\partial e_{D}}{\partial t}$ on the interface. Mathematically, we note that the effect of the Dirichlet interface error $e_{D}$ at a given time $t$ here leads to an interface error of $\frac{\partial e_{D}}{\partial t}$ in the time evolution model, which implies how the Dirichlet type of information should be properly related on different time levels when time discretization is applied, in particular when decoupling is necessary along the interface. It explains why in most of the classical explicit type of decoupled schemes, the explicit use of Dirichlet data from previous time levels, such as the Dirichlet-Dirichlet, Dirichlet-Neumann type of decoupling methods, would usually lead to stability issues. Most importantly, we note that when this scalar model problem is extended to other applications such as FSI, the Dirichlet variable will physically correspond to velocity, while the related time derivative interface terms $\rho_{1}{\bm{B}}_{1,h}\frac{\partial}{\partial t}$ in (21) will correspond to certain interface inertial force quantities due to spacial discretization, which will make the new interface condition to be derived physically meaningful in contrast to the classical Robin interface condition.

Notice that $e_{D}$ and thus $\frac{\partial e_{D}}{\partial t}$ become zero in the error equation (21) if the Dirichlet condition is enforced on the interface. Therefore, as a special case of (21), for (18) in Model $L_{h}$ we have the following interface relationship.

Theorem 3.6**.**

If $u_{h}^{{L}}$ is a solution of Model $L_{h}$, then $u_{h}^{{L}}$ satisfies the following interface condition:

[TABLE]

Mathematically, (28) is equivalent to (18) in Model $L_{h}$ since the additional two terms $\rho_{1}\left({\bm{B}}_{1,h}\frac{\partial u_{2,h}^{L}}{\partial t},\xi_{h}\right)_{\Gamma}$ and $\rho_{1}\left({\bm{B}}_{1,h}\frac{\partial u_{1,h}^{L}}{\partial t},\xi_{h}\right)_{\Gamma}$ are canceled each other due to the Dirichlet interface condition being enforced in the solution space. However, they could make significant differences numerically when further approximation is applied for time discretization as well as decoupling.

To further understand the interface mechanism implied by (28), let us assume all the regularity requirements, well-posedness, as well as error estimates for the related continuous and semi-discrete models in order to focus on illustrating the derivation and understanding the mechanism of the intrinsic Robin condition to be derived. Let us also first concentrate on the terms inside the two brackets in (28), denoted by

[TABLE]

which actually corresponds to all the terms in (18) of the lumped-mass model. Because (18) in Model $L_{h}$ is an approximation to (9) in Model $S_{h}$ due to the lumped-mass condensation of the $L^{2}$ inner product to the interface, we may relate the terms in $F$ back to (9) in Model $S_{h}$ by approximating the lumped-mass terms using $\rho_{1}\left(\frac{\partial u_{1,h}^{L}}{\partial t},\mathcal{L}_{1,h}\xi_{h}\right)_{\Omega_{1}}$ and $\rho_{2}{\left(\frac{\partial u_{2,h}^{L}}{\partial t},\mathcal{L}_{2,h}\xi_{h}\right)_{\Omega_{2}}}$ with an error of $O(h^{2})$, where the error between the inner product and its corresponding lumped-mass approximation in ${\bm{V}}_{h}$ is $O(h^{2})$ as showed in [22]. So, we have the estimate for $F$ with the solution of Model $L_{h}$, yet in terms of the expression in Model $S_{h}$ as follows.

[TABLE]

Now that Model $S_{h}$ is also the standard finite element approximation to Model C, while Model $L_{h}$ is an approximation to Model C, too. Therefore, by using (19) we may approximate $u_{1,h}^{L},~{}u_{2,h}^{L}$ in (30) by the continuous solution $u_{1}^{C},~{}u_{2}^{C}$ with the error of $O(h)$ to further relate the expression $F$ to the weak form of the coupled model, yet restricted to the test function space ${\bm{V}}_{h}$:

[TABLE]

Under sufficient regularity assumptions for the continuous solution $u^{C}$ and following the standard procedure, integrating by parts in (31) and using the strong form of the coupled PDE model leads to

[TABLE]

Therefore, we see that $F$ actually corresponds to the Newmann condition in the original interface matching conditions, or physically presenting the jump of the flux across the interface.

Combining (32) with (28), we get

[TABLE]

Approximating the continuous solution $(u_{1}^{C},~{}u_{2}^{C})$ in (33) back to the solution $(u_{1,h}^{L},~{}u_{2,h}^{L})$ of the approximate lumped-mass semi-discrete model, we find that (33) actually implies the following matching condition for $(u_{1,h}^{L},~{}u_{2,h}^{L})$ on the interface:

[TABLE]

which may be formally written as

[TABLE]

Similar to the above argument, if we condense the inner product from $\Omega_{2}$ to the interface and transmit the Dirichlet data $u_{2,h}$ to the other side of the interface, we may have the following similar interface relation

[TABLE]

First, note that in the above interface conditions (35) and (36), $\nabla u_{h}^{L}$ itself does not exist on $\Gamma$ because functions in the finite element space ${\bf V}_{h}$ have discontinuous derivatives on the interface, while the gradient of $u_{i,h}^{L}$ for $i=1,2$ is defined in the boundary layer cells in the subdomains on each side of the interface. Therefore, we denote each side of the interface in (34) by $\Gamma^{+}$ and $\Gamma^{-}$, respectively, to distinguish the interface values of a discontinuous function restricted from the corresponding interface side.

Most importantly, in contrast to the classical Robin interface condition, which simply takes a linear combination of the original Dirichlet and Neumann interface conditions, as introduced for constructing numerical methods for various decoupling purposes, this new type of interface condition instead combines the Neumann condition with the time derivative of the Dirichlet interface data, which physically describes the total exchange of heat or temperature balances on the interface and is referred to as the so-called intrinsic Robin condition. When such a scalar model problem is further extended to other applications such as FSI, as pointed out earlier, the related time derivative interface terms $\rho_{1}{\bm{B}}_{1,h}\frac{\partial}{\partial t}$ in the new type of Robin condition will correspond to certain interface inertial force quantities according to Newton’s second law, because the derivative of velocity gives acceleration, multiplying by the ”mass” factor $\rho_{1}{\bm{B}}_{1,h}$ then leads to a force. We therefore may also call this new type of interface condition generally as the inertial type Robin condition, because it is the inertial data then, instead of the velocity data, that are combined with the Neumann data from each side of the interface.

Indeed, together with the Dirichlet condition enforced in the solution space ${\bf V}_{h}$, the classical Robin condition, as well as the inertial type Robin condition (35) or (36) just introduced above, despite their approximation errors of $O(h)$, are all mathematically equivalent to the original interface condition set consisting of the Dirichlet and Neumann conditions. However, the classical Robin condition is mathematically introduced by kind of artificially combining the Dirichlet and Neumann conditions to obtain another set of mathematically equivalent interface conditions, yet lack of physical justification. Notice that, in applications like FSI, the Dirichlet data correspond to velocity, while the Neumann data correspond to interfacial force, therefore, there is a mismatch between these two different types of physical quantities. Such a physical mismatch is then kind of meredied by mathematically introducing a relaxation parameter, which is to be tuned for the purposes of numerical stability and convergence. The optimal relaxation parameter is then compensated for the mathematical effect of the physical mismatch. However, the optimal parameter analysis is often very difficult, even for certain simple model problems. The choice and tuning of the relaxation parameter with the existing numerical methods are practically very difficult and even impossible for most of the real applications, never mentioning the theoretical analysis and understanding of the numerical stability and convergence.

The inertial type of Robin condition, however, results from our error analysis in $Theorem$ 3.4 where we track how the information and its error are transmitted across the interface. The combination of the Dirichlet data, now through its time derivative instead, with the Neumann data becomes physically meaningful for applications like FSI because it characterizes how the inertial and stress forces should be balanced along the interface, which will in turn guide us in future approximation during time discretization and decoupling in order to maintain numerical stability and convergence properly.

4 Full Discretization and Decoupling

Implicit time discretization for the above semi-discrete models will lead to stable and convergent coupled methods. For illustration, let us consider the first order backward Euler scheme, although the discussions may be extended to high order schemes, such as the Crank-Nicolson Scheme. For simplicity, let us also consider linear finite elements in spacial discretization.

The standard finite element semi-discrete model is often applied to derive implicit coupled scheme by enforcing both Dirichlet and Neumann conditions on the current time level. To derive decoupled schemes, one may need to relax either Dirichlet or Neumann condition, or replace them by some other equivalent interface conditions such as the classical Robin condition or our new inertial type Robin condition. Problem $P_{h}$ illustrates how the Dirichlet data may be exchanged at certain step with the neighboring subdomain by using the computed data from the other side of the interface with the removal of the Dirichlet constraint from the solution space. For the purpose of relaxing the Neumann constraint at certain step instead of implicitly enforcing the Newmann condition as a natural interface condition imposed by Model C, one may take the local PDEs (1)-(2) in the original coupled model, multiplying by the corresponding test functions, integrating by parts, and then exchanging the Newmann information across the interface by applying the Newmann interface condition explicitly, and finally add them together to get a coupled form defined in the tensor product space of $V_{1}\times V_{2}$. This will, of course, result in the corresponding interface terms involving the Neumann data on the other side of the interface, which may later be approximated by using computed information for the purpose of decoupling. For instance, we may formally derive the following quasi-weak formulation:

**Problem EC: A Quasi-Weak Formulation of the Continuous Coupled Model

**

Find $u^{EC}=(u_{1}^{EC},u_{2}^{EC})\in{\bm{V}}$, such that

[TABLE]

Note that we call this the quasi-weak formulation because the two interface terms $\left(\beta_{2}\nabla u_{2}^{EC}\cdot{\bm{n}}_{2},v_{1}\right)_{\Gamma}$ and $\left(\beta_{1}\nabla u_{1}^{EC}\cdot{\bm{n}}_{1},v_{2}\right)_{\Gamma}$ are not defined if $u^{EC}$ is only in the $H^{1}$ space ${\bm{V}}$, unless the PDE solution has further regularity of at least $H^{3/2}$ so that the traces of $\nabla u_{1}^{EC}$ and $\nabla u_{2}^{EC}$ may be defined on the interface. However, the corresponding discrete counterparts are well defined as for the same reason as in the intrinsic or inertial Robin condition. Its corresponding finite element approximation then reads:

**Problem $ES_{h}$: A Semi-Discrete Model with Finite Element Approximation in Space

**

Find $u_{h}^{ES}=(u_{1,h}^{ES},u_{2,h}^{ES})\in{\bm{V}}_{h}$, such that

[TABLE]

For comparison, let us first review a couple of typical decoupled methods by classical decoupling approaches in the literature. Algorithm 1 describes the decoupled Dirichlet-Neumann (DN) scheme [15, 19] based on decoupling the implicit scheme from Problem $ES_{h}$ by applying the Dirichlet (3) and Neumann interface conditions (4) at alternating time levels with the use of computed data from the previous step on the other side of the interface during time marching.

In the conventional Robin-Robin approach [2, 11], the equivalent classical Robin-Robin interface conditions are introduced by combining the original Dirichlet and Neumann conditions with two mathematical relaxation parameters $\alpha_{1}$ and $\alpha_{2}$ as

[TABLE]

Algorithm 2 describes the corresponding decoupled Robin-Robin (RR) scheme by applying the two Robin conditions above at alternating time levels during time marching. As discussed earlier, such an approach is lack of physical justification, and the determination of the optimal relaxation parameters is difficult both theoretically and practically for complicated real applications.

We now propose our new decoupling approach based on the intrinsic or inertial type Robin condition. Notice that the decoupled RR Algorithm 2 corresponds to adding the two Dirichlet terms $\alpha_{1}(u_{1,h}-u_{2,h},v_{1,h})_{\Gamma}$ and $\alpha_{2}(u_{2,h}-u_{1,h},v_{2,h})_{\Gamma}$ to Model $ES_{h}$ in the classical Robin interface condition to obtain another equivalent problem and then applying the decoupling technique similar to the Dirichlet-Neumann scheme. As discussed in the previous section, the lack of physical justification of the Robin-Robin approach results in both theoretical and practical numerical difficulties. Mathematically, instead of adding the Dirichlet terms $\alpha_{1}(u_{1,h}-u_{2,h},v_{1,h})_{\Gamma}$ and $\alpha_{2}(u_{2,h}-u_{1,h},v_{2,h})_{\Gamma}$ , one may add any other term on both sides to Problem $ES_{h}$ to obtain another equivalent model without changing the solution, for instance the parameters $\alpha_{1}$ and $\alpha_{2}$ are generally arbitrary, just like in the classical fixed point iteration where there are infinite many possibilities one might try. The key is whether one may find a particular setting that leads to a better or even possibly optimal numerical algorithm. Motivated by the interface data transmission analysis in the previous section, for instance an equivalent form (28) for the lumped-mass finite element semi-discrete approximation, it is suggested that it indeed makes sense both mathematically and physically to instead add the time derivative terms $\rho_{2}\left({\bm{B}}_{2,h}(\frac{\partial u_{1,h}}{\partial t}-\frac{\partial u_{2,h}}{\partial t}),v_{1,h}\right)_{\Gamma}$ or $\rho_{1}\left({\bm{B}}_{1,h}(\frac{\partial u_{2,h}}{\partial t}-\frac{\partial u_{1,h}}{\partial t}),v_{2,h}\right)_{\Gamma}$. This leads to an equivalent semi-discrete problem of Problem $ES_{h}$:

**Problem $ESI_{h}$: An Equivalent Semi-Discrete Problem based on the Intrinsic or Inertial Robin Conditions

**Find $u_{h}^{ESI}=(u_{1,h}^{ESI},u_{2,h}^{ESI})\in{\bm{V}}_{h}$, such that

[TABLE]

Applying the backward Euler scheme gives an implicit and stable algorithm. However, this stable algorithm is coupled on the interface through the implicit approximation for the time derivative of the two interface Dirichlet data terms $\frac{\partial u_{1,h}^{n}}{\partial t}$ and $\frac{\partial u_{2,h}^{n}}{\partial t}$, as well as the Dirichlet interface condition $u_{1,h}^{n}=u_{2,h}^{n}$ enforced in the solution space.

Our framework described here will then allows one to apply various strategies to further decouple the implicit coupled scheme corresponding to Problem $ESI_{h}$, while depending on the practical physical properties in more general real applications. For instance, as motivated by the thin-wall FSI application [12], we may keep the implicit approximation for the time derivative of the interface Dirichlet data terms $\frac{\partial u_{1,h}^{n}}{\partial t}$ and $\frac{\partial u_{2,h}^{n}}{\partial t}$ on one side, say $\frac{\partial u_{1,h}^{n}}{\partial t}\approx\frac{\delta u_{1,h}^{n}}{\delta t}=\frac{u_{1,h}^{n}-{u_{1,h}^{n-1}}}{\Delta t}$, while approximate the other term by using computed data from previous time levels explicitly on the other side of the interface, say $\frac{\partial u_{2,h}^{n}}{\partial t}\approx\frac{\delta u_{2,h}^{n-1}}{\delta t}=\frac{u_{2,h}^{n-1}-{u_{2,h}^{n-2}}}{\Delta t}$. Furthermore, to decouple the Dirichlet interface condition, we may remove the enforced condition $u_{1,h}^{n}=u_{2,h}^{n}$ in the solution space, while incorporate it into the algorithm by further approximating $u_{1,h}^{n-1}$ in the implicit approximation $\frac{\delta u_{1,h}^{n}}{\delta t}=\frac{u_{1,h}^{n}-u_{1,h}^{n-1}}{\Delta t}$ by passing the Dirichlet data from the other side with $u_{2,h}^{n-1}$, which implies to explicitly enforce the Dirichlet interface condition on the previous time level instead, namely, $\frac{\delta u_{1,h}^{n}}{\delta t}=\frac{u_{1,h}^{n}-u_{1,h}^{n-1}}{\Delta t}\approx\frac{u_{1,h}^{n}-u_{2,h}^{n-1}}{\Delta t}$. This gives a decoupled intrinsic or inertial Robin-Neumann (iRN) type algorithm:

Note that it is well known from the practical experiences that the numerical instability with decoupled methods is usually caused by making use of the computed Dirichlet data explicitly from the previous time level for decoupling the Dirichlet interface condition. It is the introduction of the time derivative of the Dirichlet data in the intrinsic or inertial Robin interface condition that is not only physically justified, but also mathematically allows for the implicit use of computed Dirichlet data from the previous time level to decouple the Dirichlet condition effectively.

In addition, we have only outlined above the basic principles in our decoupling framework based on the intrinsic or inertial type Robin interface conditions. The ingredients and strategies in a particular decoupled algorithm may vary and be fine tuned, such as whether to apply one inertial or Newmann interface term and how to decouple the Dirichlet and Newmann data by explicit approximation from using computed data, on which side or both sides or on alternating time levels, which would depend on the specific mathematical and physical properties in real applications, such as the ratio of the subdomain sizes, physical parameters like $\rho,\beta$, etc. For instance, the strategy described in Algorithm 3 works particularly effectively in the thin-wall FSI application, where it is observed that the implicit interface approximation must be maintained on the bulk fluid side like $\Omega_{1}$ here, while the interface time derivation may be relaxed by applying an explicit approximation on the thin-wall structure side for decoupling. Numerical experiments will be designed and conducted in the next section to confirm the effectiveness of the proposed decoupling approach, as well as to investigate some of these questions to provide certain insight and guidance for real applications.

5 Numerical Experiments

Numerical experiments are conducted to examine the effectiveness of the proposed decoupling approach based on the intrinsic or inertial type Robin interface conditions applied to the model problem in Section 2. For illustration, let $\beta(x,y)=2+x^{2}+y^{2}$, and the source term is defined such that the exact solution of the coupled model is given by

[TABLE]

The total time is set to be $T=1s$ in all simulation.

We first examine the effectiveness of our new decoupled iRN scheme in Algorithm 3 by comparing it with the corresponding implicit coupled scheme that is known to be numerically stable and convergent. For this purpose, we simply set $\rho_{1}=\rho_{2}=1$ and $\beta_{1}=\beta_{2}=\beta$. The errors of the corresponding numerical solutions of the coupled and decoupled schemes to the exact solution are compared in Table 1 with four levels of mesh refinement, where $\Delta t=h^{2}$ in order to examine the discretization error in space. It is clearly observed that our new iRN scheme effectively decouples the computation by properly approximating the interface data using the computed solution, while still retains the same order of accuracy with the numerical solutions of the coupled scheme. The order of convergence is numerically verified because the errors are reduced by a factor of $4$ as the mesh size and time step are refined once, which indicates the second order accuracy in space and first order accuracy in time in the $L^{2}$ norm. To further examine the stability of the decoupled iRN scheme, we increase $\Delta t$ from $=O(h^{2})$ to $=O(h)$ and plot the computed solutions as well as the exact solution at $t=1$ in Figure 2, where the stability and accuracy of the decoupled iRN scheme are clearly observed. The decoupled scheme is still numerically stable as the time step size keeps increasing, although the accuracy is reduced accordingly as expected due to bigger discretization error in time as well as in approximating the interface data from previous time levels for the decoupling purpose.

In order to compare with the classical decoupled DN and RR schemes, we may also consider to apply the corresponding intrinsic or inertial Robin conditions (35) and (36) symmetrically on both sides of the interface by adding another term $\rho_{1}\left({\bm{B}}_{1,h}(\frac{\partial u_{2,h}}{\partial t}-\frac{\partial u_{1,h}}{\partial t}),v_{2,h}\right)_{\Gamma}$ to Problem $ESI_{h}$, which leads to another version of the intrinsic or inertial Robin-Robin(iRR) type of decoupled scheme:

Similar to Table 1, we compare in Table 2 the errors between the exact solution and the numerical solutions of the decoupled iRR scheme with those of the classical decoupled DN and RR schemes, as well as the implicit coupled scheme as the reference. In terms of effectiveness, it is still observed that, for simple cases where if different decoupling techniques all work, such as our new iRR scheme, as well as the classical DN scheme and RR scheme with $\alpha_{1}=10,~{}\alpha_{2}=5$, they could have comparable effective performance. However, for the classical RR method, its performance relies on, and is very sensitive to the selection of relaxation parameters as discussed earlier, which would make it difficult for theoretical analysis, and most importantly for its application and extension to complicated real problems. For instance, it is shown in the table that the accuracy gets worse as the parameters are perturbated to $\alpha_{1}=\alpha_{2}=1$. For the classical DN approach, although it works pretty well in this simple setting, it will be shown in the next experiment that it may fail, too when physical parameters vary. In the contrary, our decoupled iRR scheme is always robust and parameter-free. It is also numerically verified that the order of accuracy of the iRR scheme is similar to the iRN scheme.

Let us finally illustrate the effects of physical properties by varying the value of $\rho$ in different subdomains on different decoupling approaches of the three decoupled schemes as compared in Table 2. Two cases are tested, where Table 3 shows the results corresponding to the parameter set $\rho_{1}=10\rho_{2}$, while in Table 4 the parameter ratio is switched by $\rho_{2}=10\rho_{1}$. It is clearly seen that the classical DN approach is very sensitive to physical properties and could fail if not applied carefully, while the decoupled iRR scheme always works perfectly.

6 Concluding Remarks

We have derived an intrinsic or inertial type Robin condition for multi-modeling problems, which is justified both mathematically and physically in contrast to the classical Robin condition. Based on this new interface condition, a decoupling approach is presented for devising effective decoupled numerical methods for applications from classical parallel computing to multi-physics applications. Numerical experiments show the effectiveness and robustness of our decoupling approach, its advantages over the existing decoupling approaches, as well as the promising potential for its application to complicated real problems in science and technology. Theoretical analysis for stability and convergence is under investigation.

Acknowledgments

The authors’ research is supported in part by the Hong Kong RGC Competitive Earmarked Research Grant HKUST16301218 and NSFC (91530319,11772281). They would also like to thank Mingchao Cai and Yiyi Huang for valuable discussions.