1 Introduction
The previous papers in this series dealt with the most common differential operators such as the linear diffusion FESTUNG1 or advection FESTUNG2 as well as with different types of discontinuous Galerkin (DG) discretizations, namely the standard local DG (LDG FESTUNG1 ) or hybridized DG (HDG FESTUNG3 ). In the time span between the previous installment FESTUNG3 and the current work, several applications involving nonlinear equations were implemented using our MATLAB / GNU Octave FESTUNG FESTUNGGithub toolbox; those include the twodimensional shallowwater equations HajdukHAR2018 and mean curvature flow BungertAF2017 ; AizingerBF2018 . The present study, however, sets a much more ambitious goal: Presenting a modelcoupling interface that allows to create very complex simulation systems consisting of multiple selfcontained simulation tools and interacting by exchanging data usually in the form of solution vectors. The more specific objectives of this work include developing an abstract coupling concept capable of supporting very general types of model coupling as well as implementing and testing this concept in the framework of FESTUNG. Just as in the previous parts of this series, our performance optimized, fully vectorized implementation with a detailed documentation is freely available as an open source software.
Very few physical systems are truly isolated and thus can be modeled in a standalone fashion without the help of some more or less realistic assumptions. Often, the setting can be simplified, and action and reaction of other physical systems can be accounted for in a singlephysics model via boundary conditions, forcing terms, parametrizations, etc. However, many important problems do not lend themselves readily to such simplifications, therefore the need for coupled models. Examples of widely used multiphysics applications include, in particular, oceanatmosphereland systems constituting the staple of climate modeling or fluidstructure interaction models very common in civil engineering and medical applications. This multiphysics label can sometimes be somewhat misleading (see a discussion of this issue in a very comprehensive overview paper Keyes2013 ) and is also often used to describe interactions between different mathematical representations, or numerical schemes, or grids, or motion scales within the same ‘physics’.
For coupled multiphysics problems, one can generally identify three main classes of setups: shared domain/ multiple physics (e.g. coupled subsurface flow/geomechanics), multiple domain/shared physics (e.g. regional grid nesting in ocean or atmosphere simulations), and multiple domain/multiple physics. The first two permit a number of performancerelevant simplifications; although our new problem implementation framework accommodates either, here we focus on the third, the most general type of setup, for which a number of key aspects have to be considered:

Modeling issues: Those include the physical and mathematical consistency questions such as conservation properties, physically meaningful interface conditions, wellposedness of the mathematical problem, etc.

Numerical issues:
Different mathematical models often require specialized numerical methods; in addition, computational meshes corresponding to separate physical systems do not necessarily match at the interface (or coincide in the case of a shared domain setup). This gives rise to a number of challenges such as controlling interpolation and conservation errors between discrete solutions, matching meshes at the model interface or using, e.g., mortars
Maday1989 to bridge between subdomain meshes. Another host of rather difficult issues arises when one attempts to analyze and to improve the stability and the accuracy of the discretizations for such coupled models. These difficulties concern the stability and the convergence of spatial discretizations and – in the case of timedependent problems – can also affect the temporal convergence that may be degraded by operator splitting techniques very common in coupled applications. 
Algorithmic issues: In addition to differences in the discretization, the solution algorithms of each submodel may be such that a well thoughtout solution strategy becomes necessary for the coupled model. This aspect is particularly important for timedependent applications, where different PDE system types and resolved spatial/temporal scales may require a careful handling of the time stepping procedure. The algorithmic approaches encompass the oneway (when the data/information only flows from one model to the other) and the twoway coupling paradigms. The latter is much more general and includes three main classes of schemes:

Fully implicit (monolithic) – all subproblems together with the interface conditions are treated as one large system solved without operator splitting. This approach incurs high implementational and computational expenses but often produces the most robust and physically consistent schemes. It is the method of choice for stationary problems but it is also known to have been used for time dependent applications – specifically to deal with highly nonlinear systems of PDEs. It also tends to allow for (significantly) larger timesteps than the other coupling schemes. The main difficulties are connected with computing the Jacobi matrix for the coupled system (or solving the nonlinear system in a matrixfree fashion) – a difficult task if the submodels are complex and implemented in separate software packages.

Iterative (internal) coupling – all subsystems are solved in each coupling step alternately until a given convergence criterion is satisfied. This methodology represents a compromise between the full and noniterative couplings and is widely used in many applications. This approach allows to use a specialized software package for each submodel and provides a mechanism to control the coupling (i.e., operator splitting) error. An iterative couping requires, however, a particular care with regard to the conservation and stability properties of the discrete scheme; in addition, the operator splitting may cause convergence to a nonphysical solution if the underlying systems have, e.g., nonsmooth coefficients Keyes2013 .

Noniterative (external) coupling works similarly to the iterative coupling but forgoes iterations between the submodels. It is simple and cheap but not sufficiently robust and accurate for some applications.


Software issues: The last but not least difficulty arising when implementing coupled models concerns the software packages used for each submodel and the interface between them. The most advanced and userfriendly features are offered by generalpurpose modeling frameworks such as DUNE^{1}^{1}1https://www.duneproject.org, FEniCS^{2}^{2}2https://www.fenicsproject.org, or similar software projects. For software packages not originating from the same family, a number of specialized couplers exist that support data exchange and asynchronous execution modes, e.g. OASIS^{3}^{3}3https://portal.enes.org/oasis.
Due to their conceptual, algorithmic, and computational complexity, coupled applications have been almost exclusively a domain of large groups/companies/institutions and tend to be developed for a specific class of applications with the corresponding requirements in mind. Thus the coupled numerical models of ocean/atmosphere/land used for climate studies are usually optimized for massively parallel execution and efficient handling of large volumes of forcing and grid data; they mostly rely on a loose (external) coupling between singlephysics submodels Valcke2012 realized in a form of a dedicated piece of software such as OASIS or MCT Larson2005 . Fluidstructure interaction or groundwater flowgeomechanics models, on the other hand, are often designed from scratch in a monolithic fashion and thus solve a fully coupled nonlinear equation system for all submodels at once RinTGVT2017 .
One of the main motivations for the current work is to make the model coupling and multiphysics (in a very general sense) capability available to users without access to sophisticated software and high performance computing resources and to enable fast prototyping and testing of multiphysics applications. We limit the focus to the algorithmic and software aspects of model coupling; the numerical issues are rather straightforward here due to DG discretizations employed in all submodels; modeling issues will be presented in a brief form. As an example application to illustrate the coupling mechanism, we consider a freesurface/groundwater flow problem. Freesurface flow is represented by the 2D shallow water equations in a vertical slice (2Dv SWE – also known as primitive hydrostatic equations in a vertical slice), whereas the groundwater is modeled using Darcy’s law (also in a twodimensional vertical slice). A stability analysis of this discretization using the same interface conditions is presented in a companion paper RuppARK2018 .
The rest of this article is structured as follows: The model problem is presented in Section 2 accompanied by its LDG discretization in Section 3. The generic problem framework is the subject of Section 4, and an implementation of the model problem in the context of this framework is demonstrated in Section 5 (more implementation details are provided in A). Numerical examples can be found in Section 6, and a conclusion and outlook make up the remainder of this paper.
2 Model problem
2.1 Computational domain and mesh
Let be a finite time interval and the coupled domain consisting of a freeflow subdomain on top of the subsurface subdomain (cf. Fig. 1), both of which with compact closure and assumed to be polygonally bounded and Lipschitz. Freeflow and subsurface subdomains are separated by an interior (transition) boundary . Let be a regular family of nonoverlapping partitions of into closed elements such that . The discretization of the freeflow problem (2) requires a mesh consisting of trapezoidal elements with strictly vertical parallel edges – as depicted in Fig 1. Under this condition, the twodimensional elements are aligned in vertical columns, each corresponding to a onedimensional element on the axis later used to discretize the water height. For simplicity, the subsurface problem also uses a mesh of the same type. Denoting by the standard orthogonal projection operator onto the axis, this geometry and arrangement of mesh elements produces a nonoverlapping partition of denoted by by simply projecting elements of .
2.2 Subsurface problem
Let the boundary of be subdivided into Dirichlet and Neumann parts. We consider the timedependent Darcy equation describing water transport through fully saturated porous media, where is generally understood as the hydraulic head. The coefficient denotes the specific storativity of the porous medium and the hydraulic conductivity. Division by and setting , yields
(1a)  
where (a uniformly symmetric positive definite matrix) and are considered time and space dependent. Eq. (1a) is complemented by the following boundary and initial conditions:  
(1b)  
(1c)  
(1d) 
where denotes the outward unit normal, and is the given initial and the boundary data.
2.3 Freeflow problem
The boundary of the freeflow domain is assumed to consist of top , bottom , and lateral sections with the latter subdivided into land , open sea , river , and radiation parts. All boundaries except are regarded as timedependent, but, for brevity, the time variable is omitted. Then the 2Dv shallow water equations for velocity and total water height are given as
(2a)  
(2b)  
(2c)  
where denotes the acceleration due to gravity, is a source term, and are the freesurface elevation and the bathymetry with respect to some datum, respectively, and is a diffusion tensor satisfying the same conditions as those imposed on . The following boundary and initial conditions are specified for Eqs. (2a)–(2c):  
(2d)  
(2e)  
(2f)  
(2g)  
(2h)  
(2i) 
To summarize the different types of boundary conditions, note the following:

Bottom boundary : noslip ( with tangential vector ) and prescribed normal flow , which translates to in uncoupled simulations and interface condition (3a) in coupled simulations;

free surface and radiation boundary : vanishing normal derivative of the flow, i. e., ;

land boundary : no normal flow resulting in due to the strictly vertical boundary;

open sea boundary : prescribed water height (e. g., due to tidal forcing), i. e., and vanishing normal derivative of the flow, ;

river boundary : prescribed horizontal velocity and water height .
Note that different types of boundaries should be pieced together in a compatible manner in order to ensure wellposedness of the initialboundaryvalue problem.
2.4 Interface conditions
The coupling conditions at the interface boundary have been motivated and described in some detail in our analysis paper RuppARK2018 . They impose the continuity of the normal flux (3a) and the continuity of the dynamic pressure / head (3b).
(3a)  
(3b) 
where and denote the outward unit normals on with respect to and , correspondingly.
3 LDG Discretization
3.1 Variational formulation of the subsurface system
To formulate an LDG scheme for system (1), we first introduce an auxiliary vectorvalued unknown and rewrite (1) in mixed form also making the necessary changes to the boundary conditions:
(4a)  
(4b)  
(4c)  
(4d)  
(4e) 
The discontinuous nature of DG approximations permits formulating the variational system on an elementbyelement basis. Denoting by the unit exterior normal to , we multiply both sides of equations (4a), (4b) with smooth test functions , , correspondingly, and integrate by parts over to obtain
To improve readability, the space and time arguments are omitted – when this does not lead to ambiguities.
3.2 Variational formulation of the freeflow system
In the elevation equation (2a) and momentum equations (2b), we combine the advective fluxes in the primitive numerical fluxes
and, as in the subsurface problem, write our system in mixed form introducing an auxiliary unknown :
(5a)  
(5b)  
(5c)  
(5d)  
(5e)  
(5f)  
(5g)  
(5h)  
(5i)  
(5j) 
Similarly to Sec. 3.1, we formulate the variational system on an elementbyelement basis, multiply both sides of equations (5b)–(5d) by smooth test functions , , and , correspondingly, and integrate by parts over to obtain
For Eq. (5a), we denote by the onedimensional element corresponding to (cf. Sec. 2.3), multiply both sides by a smooth test function , and integrate by parts.
To keep the notation uniform, integrals over zerodimensional domains in denote evaluation at the respective point.
The computation of the depthintegrated velocity used in can be performed by summing the integrals over elements for which holds DawsonAizinger2005 ; AizingerPDPN2013 . Then the depthintegrated horizontal velocity on is given by
(6) 
with and being the number of elements for which . For brevity and readability we drop the “” in the following.
3.3 Definitions and preliminaries
Before describing the DG scheme for (4), (5), (6), we introduce some notation. Let , then denotes the set of interior edges, the set of boundary edges, and the set of all edges . For an interior edge shared by elements and , we define the onesided values of a scalar quantity on by
respectively. For a boundary edge , only the definition on the left is meaningful. The average and the jump of on are then given by
respectively. Note that is a vectorvalued quantity. Further notation is introduced on the first use and is summarized in the Index of notation.
For , the standard tensorproduct polynomial space of degree in each variable on , we define the broken polynomial space by
Note that we could also use instead of with some changes to the computation of the depthintegrated velocity in the free flow problem Aizinger2004 . Moreover, the projection is defined as
The corresponding componentwise generalization for vectorvalued functions is denoted by the same symbol.
3.4 Local basis representation and transformation rules
In the following, we denote by an element of and use a local numbering scheme to identify its vertices and edges , . Consequently, an interior edge belonging to neighboring elements is identified by . Due to the structure of the mesh (cf. Fig. 2), the local edge index is directly deduced from via
As in our previous works FESTUNG1 ; FESTUNG2 ; FESTUNG3 , we employ a mixture of algebraic and numerical indexing styles: for instance, means all possible combinations of a fixed element index with local edge indices such that lies in . This implicitly fixes the numerical indices which accordingly can be used to index matrices or arrays. At some points of the discretization of the freeflow problem, we require a distinction between lateral edges on the one hand and top and bottom edges on the other. For simplicity, we refer to them as ‘vertical’ and ‘horizontal’ edges and mark all sets of edges correspondingly using superscripts ‘v’ or ‘h’ although the latter are not necessarily orthogonal to the direction of gravity (cf. Fig. 2). For example, and are the sets of all vertical and horizontal interior edges, respectively. The numbering scheme of the mesh ensures that we have local edge index for horizontal and for vertical edges.
Following the standard practice, we define our discrete solution expressed in some finite basis of using a reference element. For this, we use the unit reference square as shown in Fig. 2 and specify for a diffeomorphism
We choose a set of orthonormal with respect to the inner product basis functions for , define the basis of as , and obtain representation
With explicitly defined, the mapping can be expressed in terms of the vertices of :
(7a)  
The Jacobian of the mapping and its determinant are given as  
(7b)  
(7c) 
Due to generally nonparallel top and bottom edges of the trapezoidal elements, the entries of the Jacobian matrix are not constants but rather affinelinear functions of . The local node numbering in our mesh preserves the orientation of the reference element, thus , holds. Using the Jacobian we obtain the componentwise rule for the transformation of the gradient
(8)  
For our choice of as the space of ansatz functions, we define basis functions on the reference square via tensor products of onedimensional Legendre polynomials as
(9) 
where and (see Table 1). Closedform expressions for the onedimensional basis functions on the reference interval up to order three are given by:
0  1  2  3  

1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16 
Any function implies by , i. e., , and, in particular, for all . For integrals over element and edge , we use transformation formulas
(10a)  
(10b)  
and, for integrals over onedimensional elements , we have  
(10c) 
with implying by and using the standard linear mapping
Comments
There are no comments yet.