Decoupling PDE Computation with Intrinsic or Inertial Robin Interface Condition

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.



There are no comments yet.


page 18


On the convergence and applications of the inertial-like method for null-point problems

In this paper, we propose two novel inertial-like algorithms for solving...

Solving Parabolic Moving Interface Problems with Dynamical Immersed Spaces on Unfitted Meshes: Fully Discrete Analysis

Immersed finite element (IFE) methods are a group of long-existing numer...

Wheel-Rail Interface Condition Estimation (W-RICE)

The surface roughness between the wheel and rail has a huge influence on...

An energy-based coupling approach to nonlocal interface problems

Nonlocal models provide accurate representations of physical phenomena r...

Stability Analysis of Interface Conditions for Ocean-Atmosphere Coupling

In this paper we analyze the stability of different coupling strategies ...

Numerical analysis of the Landau-Lifshitz-Gilbert equation with inertial effects

We consider the numerical approximation of the inertial Landau-Lifshitz-...

FDTD schemes for Maxwell interface problems with perfect electric conductors based on the correction function method

In this work, we propose FDTD schemes based on the correction function m...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 and be the two subdomains, with the interface as showed in Figure 1, where and

are the unit outward normal vectors to

and , respectively.

Figure 1: Computational Domains with Interface

The coupled PDE model is described as:

  in (1)
  in (2)
  on (3)
  on (4)
  on (5)
  on (6)

where and are temperature, and are density, and are material diffusivity, and 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 .

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 and , 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 , and/or , 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


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

Model C: The Weak Form of the Continuous Coupled Model

Find , such that


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 : The Semi-Discrete Model with Finite Element Approximation in Space

Find , such that


where is a standard finite element subspace of . Note now that Model 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 of Model no longer satisfies the interface conditions of the original coupled model although the Dirichlet condition remains enforced in the solution space , 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 be a matrix and be the corresponding lumped-mass matrix




Denote consists of continuous, piecewise linear functions on a quasiuniform family of triangulations of any domain , where or . The inner product in is defined by

Let be a triangle of the triangulation and be its vertices. Define the quadrature formula for any function


then define the lumped-mass approximation of inner product by


Let be the finite element basis, the following two properties hold [22]:

  • for ,

  • .

From the properties, if is generated by , then is generated by .

Denote be the corresponding discrete spaces with continuous piecewise linear functions. Applying the lumped-mass technique to , we may build connections among the inner products in and . Define two discrete lifting operators


such that the nodal values of and vanish out of for all . The interface operator is defined by , where is the adjoint operator of w.r.t the lumped-mass inner product in . is also defined similarly. Based on the discrete operators introduced above, we have


The following Lemmas are useful in the analysis below. For simplicity, we denote by .

Lemma 1

For all , we have



For , it can be decomposed into with . From the definition of lumped-mass inner product (13), we can conclude that for all since all the nodal values of are . Therefore, using (15), we have

Similarly, we also have

Lemma 2

For all , we have


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

Model : The Lumped-Mass Semi-Discrete Model

Find , such that


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 has regularity, for each subdomain we have


There are many important applications in numerical computation where the Dirichlet condition enforced in the solution space 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 by removing the Dirichlet interface condition from the solution space , where the modification of is defined as .

Problem : A Perturbation to Model

Consider a function satisfying


The following theorem is revealing.

Theorem 1

If satisfies (20) in Problem , then satisfies the following condition:



Given , taking and in (20), we have


The lumped-mass property (16) in Lemma 3.1 allows us to condense the inner product in subdomain to the interface :


In order to pass the Dirichlet information from to , denote , we have


Applying (24) to (22) gives


Reorganizing the terms by separating the information on each subdomain yields


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 for subdomain , and the time derivative term being condensed to the interface from subdomain , but then exchanged the Dirichlet interface data with subdomain 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 to equation (26), which gives


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 1 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 on the interface. Mathematically, we note that the effect of the Dirichlet interface error at a given time here leads to an interface error of 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 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 and thus 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 we have the following interface relationship.

Theorem 2

If is a solution of Model , then satisfies the following interface condition:


Mathematically, (28) is equivalent to (18) in Model since the additional two terms and 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


which actually corresponds to all the terms in (18) of the lumped-mass model. Because (18) in Model is an approximation to (9) in Model due to the lumped-mass condensation of the inner product to the interface, we may relate the terms in back to (9) in Model by approximating the lumped-mass terms using and with an error of , where the error between the inner product and its corresponding lumped-mass approximation in is as showed in [22]. So, we have the estimate for with the solution of Model , yet in terms of the expression in Model as follows.


Now that Model is also the standard finite element approximation to Model C, while Model is an approximation to Model C, too. Therefore, by using (19) we may approximate in (30) by the continuous solution with the error of to further relate the expression to the weak form of the coupled model, yet restricted to the test function space :


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


Therefore, we see that 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


Approximating the continuous solution in (33) back to the solution of the approximate lumped-mass semi-discrete model, we find that (33) actually implies the following matching condition for on the interface: