Advances in Mathematical Physics

Volume 2015 (2015), Article ID 787198, 13 pages

http://dx.doi.org/10.1155/2015/787198

## The Interaction of Iteration Error and Stability for Linear Partial Differential Equations Coupled through an Interface

^{1}Department of Mathematics, Colorado State University, Fort Collins, CO 80523, USA^{2}Department of Statistics, Colorado State University, Fort Collins, CO 80523, USA^{3}Tech X Corporation, Boulder, CO 80303, USA

Received 21 October 2014; Revised 16 December 2014; Accepted 18 December 2014

Academic Editor: Soheil Salahshour

Copyright © 2015 B. Sheehan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We investigate properties of algorithms that are used to solve coupled evolutionary partial differential equations posed on neighboring, nonoverlapping domains, where the solutions are coupled by continuity of state and normal flux through a shared boundary. The algorithms considered are based on the widely used approach of iteratively exchanging boundary condition data on the shared boundary at each time step. There exists a significant and sophisticated numerical analysis of such methods. However, computations for practical applications are often carried out under conditions under which it is unclear if rigorous results apply while relatively few iterations are used per time step. To examine this situation, we derive exact matrix expressions for the propagation of the error due to incomplete iteration that can be readily evaluated for specific discretization parameters. Using the formulas, we show that the universal validity of several tenants of the practitioner’s conventional wisdom are not universally valid.

#### 1. Introduction

The class of multiphysics problems in which one physical process operates in one domain while a second physical process operates in a neighboring domain, and the solutions in the component subdomains are coupled by continuity of state and normal flux through a common boundary, which is central to a number of applications [1]. Examples include conjugate heat transfer between a fluid and solid [2–5], Stokes-Darcy flow in karst aquifers [6] and in catalytic converters [7], modeling of the core and edge plasma flows in a tokamak reactor [8, 9], and flow in heterogeneous porous media [10–13].

In applications, the processes in the component subdomains are often themselves represented by complex, multiscale, and multiphysics models and the component models are solved with different discretization methods utilizing significantly different scales. Such coupled problems present enormous challenges for efficient high performance computing [1]. There are very strong incentives to use existing “legacy” codes for the component models and to treat component physics solvers as “black boxes.” For these reasons the most commonly encountered solution strategy is a simple iterative approach [1–4, 14–17]. In this approach, the models in each of the component subdomains are associated with one of the two boundary coupling conditions and subsequently solved independently except for data from the coupling through the boundary. The coupling data for the boundary conditions for one component model are provided by the solution of the other component model from the previous iteration. In a time dependent problem, this iteration is carried out on each time step.

This type of coupled physics problem has been extensively studied by the numerical analysis community, with the goals of deriving accurate, stable numerical methods, efficient, and accurate coupling algorithms and deriving rigorous a priori analyses, for example, in [10–15, 17, 18], as well as a posteriori analyses [2–4, 19]. However, in practice, as mentioned, the component physics is usually complex and there are large differences in the discretization strategies employed in the component subdomains necessitating significant processing of transferred information. Thus, the application of current rigorous mathematical methods and analysis is problematic if not impossible. This is partially the reason that there is a very signicant gulf between the kinds of practices carried out in high performance computational science and the “best practices” determined by the mathematical community [3].

In addition, the computational strategies employed in practice are often rationalized by a body of conventional wisdom which asserts that (1) stability is equivalent to accuracy and (2) the use of unconditionally stable implicit solvers for the component physics generally stabilizes the entire coupled simulation as long as the component solutions are stable. This is often rationalized by some sort of (Neumann) linear stability analysis applied to the simplied case of two coupled linear parabolic problems. Based on these ideas, simulations are deemed acceptable provided they do not exhibit obvious instabilities such as unbounded behavior or rapid oscillation. Moreover, due to the prohibitive computational cost, such conclusions are often based on the single computation at hand, rather than through an exhaustive study.

This paper attempts to address the difficulties that face extending rigorous convergence analysis to the complex models and discretizations that occur in practice and illustrate issues arising from the conventional wisdom. Instead of focusing on deriving conditions that guarantee good behavior as is usual for a standard convergence analysis, we adopt a different approach and develop a rigorous computable error representation that can be evaluated for any given choice of discretization parameters, thus allowing the conventional wisdom to be verified or not. We consider the canonical problem of two linear coupled parabolic equations and formulate the iterative solution method as a fixed point iteration for a block system on the discrete level. We then derive an exact formula for the error in the iterated solution. This formula is related to the Neumann series for the inverse of the full system matrix. The argument has several virtues: it is elementary, it subsumes any sort of ad hoc linear stability analysis, and it is general in the sense that it holds for a variety of discretization methods and range of scales for the component physics. Importantly, it allows for easy evaluation of various discretization parameter choices, for example, step size, space mesh, and number of iterations per time step.

We then present a detailed study of the canonical problem of two linear coupled parabolic equations that amply demonstrates shortcomings of the conventional wisdom. Firstly, we are able to dispel the notion that stability implies accuracy. In particular, we demonstrate that a divergent iteration can be part of a stable algorithm if the number of iterations used is low, while the resulting approximation is inaccurate. In this case, the user might not be aware that the iteration is divergent if they only consider whether or not obvious instability occurs. This case is particularly interesting since a serious problem is being masked by seemingly “reasonable” results. Secondly, we demonstrate that there is no guarantee that an “unconditionally stable” time integration scheme like backward Euler remains unconditionally stable if the system is not solved exactly at each time step. We explore cases where using a limited number of iterations leads to a conditional stability, which depends on the size of the time step. The values of time steps which provide stability can occupy a range with a maximum and minimum value, rather than just a maximum.

The remainder of this paper is organized as follows. Section 2 introduces the problem and the notation associated with discretization and iteration. Section 3 derives the primary results of the paper regarding the stability of the iterative solution. Section 4 provides numerical examples. Section 5 addresses the multirate case, in which the time step in one component is an integer multiple of the time step in the other component. Multirate numerical examples are included. A brief conclusion is given in Section 6.

#### 2. Problem Formulation and Discretization

We consider a system posed on a domain consisting of two neighboring, nonoverlapping, convex, and polygonal subdomains and in 1, 2, or 3 spatial dimensions that share a common polygonal interface boundary (see [1, 2, 8, 9] for applications). We illustrate in Figure 1.