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 non-linear equations were implemented using our MATLAB / GNU Octave FESTUNG FESTUNGGithub toolbox; those include the two-dimensional shallow-water equations HajdukHAR2018 and mean curvature flow BungertAF2017 ; AizingerBF2018 . The present study, however, sets a much more ambitious goal: Presenting a model-coupling interface that allows to create very complex simulation systems consisting of multiple self-contained 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 single-physics 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 multi-physics applications include, in particular, ocean-atmosphere-land systems constituting the staple of climate modeling or fluid-structure interaction models very common in civil engineering and medical applications. This multi-physics 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 multi-physics 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 performance-relevant 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, well-posedness of the mathematical problem, etc.
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., mortarsMaday1989 to bridge between sub-domain 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 time-dependent 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 sub-model may be such that a well thought-out solution strategy becomes necessary for the coupled model. This aspect is particularly important for time-dependent 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 one-way (when the data/information only flows from one model to the other) and the two-way coupling paradigms. The latter is much more general and includes three main classes of schemes:
Fully implicit (monolithic) – all sub-problems 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 non-linear systems of PDEs. It also tends to allow for (significantly) larger time-steps than the other coupling schemes. The main difficulties are connected with computing the Jacobi matrix for the coupled system (or solving the non-linear system in a matrix-free fashion) – a difficult task if the sub-models are complex and implemented in separate software packages.
Iterative (internal) coupling – all sub-systems are solved in each coupling step alternately until a given convergence criterion is satisfied. This methodology represents a compromise between the full and non-iterative couplings and is widely used in many applications. This approach allows to use a specialized software package for each sub-model 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 non-physical solution if the underlying systems have, e.g., non-smooth coefficients Keyes2013 .
Non-iterative (external) coupling works similarly to the iterative coupling but forgoes iterations between the sub-models. 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 sub-model and the interface between them. The most advanced and user-friendly features are offered by general-purpose modeling frameworks such as DUNE111https://www.dune-project.org, FEniCS222https://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. OASIS333https://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 single-physics sub-models Valcke2012 realized in a form of a dedicated piece of software such as OASIS or MCT Larson2005 . Fluid-structure interaction or groundwater flow-geomechanics models, on the other hand, are often designed from scratch in a monolithic fashion and thus solve a fully coupled non-linear equation system for all sub-models at once RinTGVT2017 .
One of the main motivations for the current work is to make the model coupling and multi-physics (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 multi-physics 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 sub-models; modeling issues will be presented in a brief form. As an example application to illustrate the coupling mechanism, we consider a free-surface/groundwater flow problem. Free-surface 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 two-dimensional 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 free-flow subdomain on top of the subsurface subdomain (cf. Fig. 1), both of which with compact closure and assumed to be polygonally bounded and Lipschitz. Free-flow and subsurface subdomains are separated by an interior (transition) boundary . Let be a regular family of non-overlapping partitions of into closed elements such that . The discretization of the free-flow problem (2) requires a mesh consisting of trapezoidal elements with strictly vertical parallel edges – as depicted in Fig 1. Under this condition, the two-dimensional elements are aligned in vertical columns, each corresponding to a one-dimensional 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 non-overlapping 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 time-dependent 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
|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:|
where denotes the outward unit normal, and is the given initial and the boundary data.
2.3 Free-flow problem
The boundary of the free-flow 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 time-dependent, but, for brevity, the time variable is omitted. Then the 2Dv shallow water equations for velocity and total water height are given as
|where denotes the acceleration due to gravity, is a source term, and are the free-surface 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):
To summarize the different types of boundary conditions, note the following:
Bottom boundary : no-slip ( 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 well-posedness of the initial-boundary-value 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).
where and denote the outward unit normals on with respect to and , correspondingly.
3 LDG Discretization
3.1 Variational formulation of the subsurface system
The discontinuous nature of DG approximations permits formulating the variational system on an element-by-element 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 free-flow system
and, as in the subsurface problem, write our system in mixed form introducing an auxiliary unknown :
Similarly to Sec. 3.1, we formulate the variational system on an element-by-element basis, multiply both sides of equations (5b)–(5d) by smooth test functions , , and , correspondingly, and integrate by parts over to obtain
To keep the notation uniform, integrals over zero-dimensional domains in denote evaluation at the respective point.
The computation of the depth-integrated velocity used in can be performed by summing the integrals over elements for which holds DawsonAizinger2005 ; AizingerPDPN2013 . Then the depth-integrated horizontal velocity on is given by
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 one-sided 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 vector-valued quantity. Further notation is introduced on the first use and is summarized in the Index of notation.
For , the standard tensor-product 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 depth-integrated velocity in the free flow problem Aizinger2004 . Moreover, the -projection is defined as
The corresponding component-wise generalization for vector-valued 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 free-flow 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 :
|The Jacobian of the mapping and its determinant are given as|
Due to generally non-parallel top and bottom edges of the trapezoidal elements, the entries of the Jacobian matrix are not constants but rather affine-linear functions of . The local node numbering in our mesh preserves the orientation of the reference element, thus , holds. Using the Jacobian we obtain the component-wise rule for the transformation of the gradient
For our choice of as the space of ansatz functions, we define basis functions on the reference square via tensor products of one-dimensional Legendre polynomials as
where and (see Table 1). Closed-form expressions for the one-dimensional basis functions on the reference interval up to order three are given by:
Any function implies by , i. e., , and, in particular, for all . For integrals over element and edge , we use transformation formulas
|and, for integrals over one-dimensional elements , we have|
with implying by and using the standard linear mapping