Time-dependent Dirichlet conditions in finite element discretizations

Abstract

For the modeling and numerical approximation of problems with time-dependent Dirichlet boundary conditions, one can call on several consistent and inconsistent approaches. We show that spatially discretized boundary control problems can be brought into a standard state space form accessible for standard optimization and model reduction techniques. We discuss several methods that base on standard finite element discretizations, propose a newly developed problem formulation, and investigate their performance in numerical examples. We illustrate that penalty schemes require a wise choice of the penalization parameters in particular for iterative solves of the algebraic equations. Incidentally, we confirm that standard finite element discretizations of higher order may not achieve the optimal order of convergence in the treatment of boundary forcing problems and that convergence estimates by the common method of manufactured solutions can be misleading.

Main article text

INTRODUCTION

In practical applications, see [1,2] for examples on flow control, a system is typically controlled via actuations at an interface. The mathematical model to use is, thus, a partial differential equation (PDE) with respect to space and possibly time posed on a domain and controls acting at the boundary. Depending on the application, the control may appear as a Dirichlet or a Neumann or Robin boundary condition.

Despite their importance in the modeling of control setups, cf. [3, Ch. 1], time-dependent inhomogeneous Dirichlet conditions have sparsely been investigated in terms of analysis and numerical approximation. Also for the elliptic or time-independent case, in textbooks on optimal control of PDEs, inhomogeneous Dirichlet conditions are often not considered because they are not of variational type, i.e., the equations are not posed in a dual space of the solution space, see, e.g., Refs. [4, Ch. 2] and [5, Ch. 2.3]. Another rather obvious obstacle is that a standard choice of trial and test function formulations implies a certain smoothness of the boundary data which may be impractical [5, Ch. 2.3].

For a general overview of the functional analysis for parabolic systems with Dirichlet boundary control, we refer to Refs. [4,6]. One basic approach is to transpose the involved elliptic operator so that the boundary conditions appear in the dynamic equations. This approach considers test functions of higher regularity and allows for rough solutions and boundary values. In the books mentioned, this method is referred to as Method of Transposition.

More recently, in the literature on numerical approximation of this type of solutions, the term very or ultra weak solutions has been used. The elliptic case is treated in Refs. [7–9], and time-dependent formulations are considered in Refs. [10,11]. The existence and the approximation of very weak solutions are well understood [7]. The key ingredient is the proper approximation of functions at the boundary via a projection [7,8,10].

An alternative approach of relaxing the boundary constraint via a penalization term in Robin boundary conditions has been investigated in Refs. [12,13].

The scope of the work presented is the assessment of the numerical treatment of boundary control problems in view of employing standard finite dimensional state space system theory for optimal control and model reduction; see Ref. [14] for an application example. The main criterion is that we can use standard continuous Galerkin schemes and that the spatially discretized problem can be written in the form:

(1)$$\dot{v}(t)=g(t,v,u)$$

or, in the linear case, in the form:

(2)$$\dot{v}(t)=A(t)v(t)+B(t)u(t).$$

We will consider algebraic manipulations of spatial discretizations of the standard formulation, as well as reformulations of the abstract equations and discuss their performance in numerical approximation of convection-diffusion problems. Apart from the value of an overview and a comparison of more or less well-known approaches, this paper provides evidence and insight into two phenomena that are important for the numerical analysis but that have not gained particular attention yet:

1. The convenient and analytically well-understood approach of the approximative Robin boundary conditions will likely fail if the state equations are solved only up to a given relative residual.

2. In the considered example, the convergence order of standard finite element schemes of polynomial degree 2 for time-dependent boundary-driven problems is lower than what one would expect from the convergence order for stationary problems. This lower convergence rate is not detected by the method of manufactured solution that is often used to numerically determine the convergence.

In this manuscript, we define consistency, i.e., the reformulation does not change the solution, on the semi-discrete level. Hence, we take the point of view that the solution of the equivalent representation will converge, if the chosen discretization scheme converges. However, this might not be the case, see Ref. [15, Ch. 1] for an example considering the Navier–Stokes equations. In short, the consistency of the algebraic manipulations with reformulations of the abstract equations is of highest importance for stable and convergent approximations. We will consider this issue for the treatment of Dirichlet conditions separately in a forthcoming paper.

This paper is organized as follows. In the section Generic problem formulations, we state the type of problems that we will consider both in an abstract setting and after a spatial discretization. In the section Rewriting the spatially discretized equations, we consider approaches that reformulate the spatially discretized equations into the desired form. In the section Incorporation via variational formulations and their discretizations, we discuss methods that reformulate the abstract equations such that a spatial discretization is a system of distributed type (1). In the section Numerical tests, we report on numerical tests concerning the approximation properties of the introduced methods. We conclude the paper by summarizing remarks and an outlook.

GENERIC PROBLEM FORMULATIONS

We will define a general continuous formulation that covers weak formulations of many PDEs from the modeling of physical phenomena. Also, we state the generic form of a spatial semi-discretization. We will restrict the considerations to the scalar case.

Continuous equations

Let $\Omega \in {ℝ}^d,d\in \{2,3\}$, be a bounded and regular domain such that the trace theorem as formulated in Ref. [16, Thm. 3.1] applies. Let Γ be its boundary. We define the Sobolev spaces $\mathcal{V}:=W^{1,2}(\Omega )$ and $\mathscr{H}:=L^2(\Omega )$ and the dual space $\mathcal{V}\prime$ of $\mathcal{V}$ with respect to the continuous embedding of $\mathcal{V}$ in $\mathscr{H}$ to get:

$$\mathcal{V}\hookrightarrow \mathscr{H}\hookrightarrow \mathcal{V}\prime .$$

We also introduce abbreviations for the trace spaces, cf. [17, Ch. 1.1], via

$$\mathcal{Q}\prime =W^{\frac{1}{2},2}(\Gamma )\text{and}\mathcal{Q}=\mathcal{Q}\prime \prime :=ℒ(\mathcal{Q}\prime ,ℝ),$$

the space of bounded linear functionals on $\mathcal{Q}\prime$.

Let

(3)$$\gamma :\mathcal{V}\to \mathcal{Q}\prime$$

be the trace operator as defined, e.g., in Ref. [17, Thm. 1.5].

We state the prototype of the continuous problem.

Problem 2.1. Let T > 0 and consider $\mathcal{A}:(0,T)\times \mathcal{V}\to \mathcal{V}\prime$. For $ℱ\in L^2(0,T;\mathcal{V}\prime )$, for ${\upsilon }_0\in \mathscr{H}$, and $\mathcal{U}\in L^2(0,T;\mathcal{Q}\prime )$, find $\upsilon$ with $\upsilon (t)\in \mathcal{V}$ and $\dot{\upsilon }(t)\in \mathcal{V}\prime$, a.e. on (0, T), so that:

(4a)$$\dot{\upsilon }(t)-\mathcal{A}(t,\upsilon (t))=ℱ(t),$$

(4b)$$\gamma \upsilon (t)=\mathcal{U}(t),$$

holds for almost all $t\in (0,T)$, and so that $\upsilon (0)={\upsilon }_0$ in $\mathscr{H}$.

The system of abstract Equations (1) contains common weak formulations of PDEs that model physical phenomena, cf. [18]. We will not address time regularity here and, thus, leave the properties of the mappings $t\mapsto \mathcal{A}(t,\upsilon (t))$ and, e.g., $t\mapsto \dot{\upsilon }(t)$ undefined in the statement of Problem 2.1.

As an example, we consider the convection diffusion equation that models the propagation of a scalar quantity ρ due to convection and diffusion in a domain.

Problem 2.2. Given a domain $\Omega \in {ℝ}^d$, a diffusion parameter $\nu$, a convection wind β, with $\beta \left(x,t\right)\in {ℝ}^d$ for time t > 0 and $x\in \Omega$, an initial value ρ₀, and a function $g$, with $g\left(t\right):\Gamma \to ℝ$ prescribing the boundary conditions, find a function ρ of space and time that satisfies:

(5a)$$\dot{\rho }(t)+\beta ·\nabla \rho (t)-\nu \Delta \rho (t)=0,$$

(5b)$$\rho {|}_{\Gamma }(t)=g(t),$$

and ρ(0) = ρ₀.

In standard weak formulations, assuming $\upsilon \in \mathcal{V}:=W^{1,2}(\Omega )$, Problem 2.2 is of the type of Problem 2.1, with, e.g., $\mathcal{A}$ defined via:

(6)$${\langle \mathcal{A}(t,\upsilon (t)),\varphi \rangle }_{\mathcal{V}\prime ,\mathcal{V}}={\displaystyle {\int }_{\Omega }(w·\nabla \upsilon (t),\varphi )}+\nu (\nabla \upsilon (t),\nabla \varphi )\text{d}\omega -{\displaystyle {\int }_{\Gamma }\nu }\left(\frac{\partial \upsilon }{\partial n}(t),\varphi \right)\text{d}\gamma ,$$

for all $\varphi \in \mathcal{V}$ and with $\partial /\partial n$ denoting the normal derivative. Here and in what follows, the pairing $(·,·)$ denotes the inner product in the spaces under consideration. The boundary integral in Equation (6) is only well defined for sufficiently regular solutions and test functions. For $\upsilon (t)$, $\varphi \in V=W^{1,2}(\Omega )$, it holds that $\frac{\partial \upsilon }{\partial n}(t){|}_{\Gamma }\in W^{-\frac{1}{2},2}(\Gamma )$ and $\varphi {|}_{\Gamma }\in W^{\frac{1}{2},2}(\Gamma )$, so that the term ${\displaystyle {\int }_{\Gamma }\nu }\left(\frac{\partial \upsilon }{\partial n}(t),\varphi \right)\text{d}\gamma$ is well defined as the continuous extension of the inner product in L²(Γ), see Ref. [17, Ch. 1.1].

Note that there are other possible choices for a weak formulation [11].

The boundary condition in Equation (4), viewed as a constraint, can also be incorporated using the dual operator of $\gamma :\mathcal{V}\to \mathcal{Q}\prime$ and a so-called Lagrange multiplier. Then, under certain smoothness and consistency conditions [19], Problem 2.1 is equivalent to:

Problem 2.3. Let T > 0 and consider $\mathcal{A}:(0,T)\times \mathcal{V}\to \mathcal{V}\prime$. For $ℱ\in L^2(0,T;\mathcal{V}\prime )$, ${\upsilon }_0\in \mathscr{H}$, and $\mathcal{U}\in L^2(0,T;\mathcal{Q}\prime )$, find $\upsilon$ with $\upsilon (t)\in \mathcal{V}$ and $\dot{\upsilon }(t)\in \mathcal{V}\prime$ and Λ with $\Lambda (t)\in \mathcal{Q}$, a.e. on (0, T), so that:

(7a)$$\dot{\upsilon }(t)-\mathcal{A}(t,\upsilon (t))-\gamma \prime \Lambda (t)=ℱ(t),$$

(7b)$$\gamma \upsilon (t)=\mathcal{U}(t),$$

hold for almost all $t\in (0,T)$, and so that $\upsilon (0)={\upsilon }_0$ in $\mathscr{H}$.

Note that, in general, the Lagrangian multiplier resides in the dual space of the constraint. In the considered case, where $Q\prime =W^{\frac{1}{2},2}(\Gamma )$ is a Hilbert space, we can deliberately identify $Q\prime \prime =Q$.

Spatially discretized equations

We consider a generic spatial discretization of the introduced equations. Let $V\subset \mathcal{V}$ be a finite dimensional subspace spanned by the basis functions ${\left\{{\psi }_i\right\}}_{i=1}^{n_v}$. As it is standard for spatial discretizations of PDEs, we consider nodal bases, i.e., the basis functions are associated with nodes of a mesh and they have local support. We consider the decomposition:

$$V=V_I\oplus V_{\Gamma },$$

where $V_I=\text{span}{\left\{{\psi }_i\right\}}_{i=1}^{n_I}$ is the space spanned by the basis functions that are associated with nodes in the inner and that are zero at the boundary. Accordingly, $n_I$ is the number of nodes in the inner and $V_{\Gamma }\subset V$ is spanned by the basis functions ${\left\{{\psi }_i\right\}}_{i=n_I+1}^{n_v}$ that have nonzero values at the boundary. We will use the abbreviation dof to address a degree of freedom that is represented by a basis function of V. Note that the considered splitting of V is not necessarily orthogonal.

Thus, at time t, the function $\upsilon (t)\in \mathcal{V}$ is to be approximated by a finite dimensional function $\upsilon (t)\in V$ or the vector $\upsilon (t)\in {ℝ}^{n_v}$ containing the coefficients of the expansion in the considered basis. We will assume that $\upsilon =({\upsilon }_I,{\upsilon }_{\Gamma })$ is partitioned, with ${\upsilon }_I$ being associated with V_I and ${\upsilon }_{\Gamma }$ being associated with V_Γ, i.e., the parts of V that live in the inner and at the boundary of the considered domain.

Without further mentioning, for a function $\upsilon \in V$, we will identify ${\upsilon }_I$ and ${\upsilon }_{\Gamma }$ with their coefficient vectors of the expansion in Equation (8) and simply write:

(8)$$\upsilon ={\upsilon }_I+{\upsilon }_{\Gamma }={\displaystyle \sum _{i=1}^{n_I}{\upsilon }_i}{\psi }_i+{\displaystyle \sum _{i=n_I+1}^{n_{\upsilon }}{\upsilon }_i}{\psi }_i.$$

We will consider test spaces that are subspaces of V. If only Dirichlet conditions are posed, the standard test space is V_I. Otherwise, all boundary dofs that are not fixed by a Dirichlet condition are included in the test space.

Generally, in the assembled coefficient matrices, rows will correspond to dofs in the test space and columns will correspond to dofs in the ansatz space. In particular, we will consider complying partitions of the coefficient matrices like the mass matrix:

$$M:={\left[\begin{array}{c}{({\psi }_i,{\psi }_j)}_{\mathscr{H}}\end{array}\right]}_{i,j=1,\dots ,n_v}$$

with respect to the test space,

$$M=\left[\begin{array}{c}{M}_I\\ M_{\Gamma }\end{array}\right],$$

and, once more, with respect to the trial space,

(9)$$M_I=\left[\begin{array}{cc}{M}_{II}& M_{I\Gamma }\end{array}\right],$$

where

$$M_{II}:={[{({\psi }_i,{\psi }_j)}_{\mathscr{H}}]}_{i,j=1,\dots ,n_I}\text{and}{M}_{I\Gamma }:={[{({\psi }_i,{\psi }_j)}_{\mathscr{H}}]}_{i=1,\dots ,n_I,j=n_I+1,\dots ,n_{\upsilon }}$$

are the parts associated with the inner dofs and the part of the mass matrix that relates to the boundary dofs tested against the inner nodes, respectively.

Similarly, we define the discrete approximation $A:(0,T)\times {ℝ}^{n_{\upsilon }}\to {ℝ}^{n_{\upsilon }}$ to $\mathcal{A}$ as:

$$A(t,\upsilon )={\left[\begin{array}{c}{\langle \mathcal{A}(t,\upsilon ),{\psi }_i\rangle }_{\mathcal{V}\prime ,\mathcal{V}}\end{array}\right]}_{i=1,\dots ,n_{\upsilon }}$$

and $A_I:(0,T)\times {ℝ}^{n_{\upsilon }}\to {ℝ}^{n_I}$ as its restriction to the test functions of the inner nodes, where, again, we have associated a vector $\upsilon \in {ℝ}^{n_{\upsilon }}$ with a function in V via (8). If $\mathcal{A}$ is linear, then its approximation on $A:(0,T)\times {ℝ}^{n_{\upsilon }}\to {ℝ}^{n_{\upsilon }}$ can be assembled as a matrix-valued function via:

$$A(t)={\left[\begin{array}{c}{\langle \mathcal{A}(t,{\psi }_j),{\psi }_i\rangle }_{\mathcal{V}\prime ,\mathcal{V}}\end{array}\right]}_{i,j=1,\dots ,n_{\upsilon }},$$

with the partitions A_I, A_II, and A_IΓ as they were defined for _M in Equation (9).

The discrete approximation $f:(0,T)\to {ℝ}^{n_{\upsilon }}$ to the right-hand side $ℱ:(0,T)\to \mathcal{V}\prime$ is given as:

$$f(t)={\left[\begin{array}{c}{\langle ℱ(t),{\psi }_i\rangle }_{\mathcal{V}\prime ,\mathcal{V}}\end{array}\right]}_{i=1,\dots ,n_{\upsilon }}.$$

We will not distinguish notationally between f and its restriction to the inner test functions.

To assign the boundary values, we simply assign the dofs associated with the corresponding boundaries via:

(10)$$G\upsilon =u,\text{where }G=\left[\begin{array}{cc}0& I\end{array}\right]\in {ℝ}^{n_{\upsilon }-n_I,n_{\upsilon }}$$

and where $u\in {ℝ}^{n_{\upsilon }-n_I}$ is a vector that will contain the current value of the boundary control at the given locations. As defined in Equation (10), the operator $G$ picks out the boundary dofs of a function $\upsilon \in V$ and assigns the control values. Note, however, that there are other discrete approximations to the trace operator $\gamma$ that consider, e.g., the inner product in $Q\prime$ or include suitable projections [7,8].

Thus, if we assume $\upsilon (t)\in V$ and if we test against the basis functions of V_I, the generic spatial discretization of Problem 2.1, that treats the boundary separately from the differential equation is of the form:

Problem 2.4. Let T > 0, $n_{\upsilon }$, $n_I\in \mathbb{N}$, and $n_d:=n_{\upsilon }-n_I$ and consider $A_I:(0,T)\times {ℝ}^{n_{\upsilon }}\to {ℝ}^{n_I}$, $G\in {ℝ}^{n_d,n_{\upsilon }}$, and $M_I\in {ℝ}^{n_I,n_v}$ as defined in the beginning of the section, Spatially discretized equations. For $f\in L^2(0,T;{ℝ}^{n_I})$, $\alpha \in {ℝ}^{n_{\upsilon }}$, and $u\in L^2(0,T;{ℝ}^{n_d})$ find $\upsilon$ with $\upsilon (t):(0,T)\to {ℝ}^{n_{\upsilon }}$, so that:

(11a)$$M_I\dot{\upsilon }(t)-A_I(t,\upsilon (t))=f(t),$$

(11b)$$G\upsilon (t)=u(t),$$

hold for almost all $t\in (0,T)$ and $\upsilon (0)=\alpha$.

For the system of Problem 2.3 with the multiplier, a possible spatial discretization defines a differential equation considering also the boundary parts, cf. [19,20]. It generically takes the form:

Problem 2.5. Let $T>0$, $n_{\upsilon }$, $n_I\in \mathbb{N}$, and $n_d:=n_{\upsilon }-n_I$ and consider $A:(0,T)\times {ℝ}^{n_{\upsilon }}\to {ℝ}^{n_{\upsilon }}$, $G\in {ℝ}^{n_d,n_{\upsilon }}$, and $M\in {ℝ}^{n_{\upsilon },n_{\upsilon }}$ as defined in the beginning of the section Spatially discretized equations. For $f\in L^2(0,T;{ℝ}^{n_{\upsilon }})$, $\alpha \in {ℝ}^{n_{\upsilon }}$, and $u\in L^2(0,T;{ℝ}^{n_d})$ find $\upsilon :(0,T)\to {ℝ}^{n_{\upsilon }}$ and $\lambda :(0,T)\to {ℝ}^{n_d}$, so that:

(12a)$$M\dot{\upsilon }(t)-A(t,\upsilon (t))-G^{\text{T}}\lambda (t)=f(t),$$

(12b)$$G\upsilon (t)=u(t),$$

hold for almost all $t\in (0,T)$ and $\upsilon (0)=\alpha$.

For illustration purposes, we will use the linear time-invariant case of Problem 2.4, for which A_I is a linear map given as a matrix $A_I\in {ℝ}^{n_I,n_{\upsilon }}$ and write (4) as:

(13a)$$M_I\dot{\upsilon }(t)-A_I\upsilon (t)=f(t)$$

(13b)$$Gv(t)=u(t).$$

More often than not, we will omit the time dependency of the variables and functions.

Remark 2.6. Until now we have not addressed time regularity, but, for sufficiently smooth input functions, we expect to obtain solutions in the classical sense. Only the values at the boundaries may have a jump at t = 0, since consistency with the boundary conditions is not possible for an arbitrary input. This is in line with the infinite dimensional setting, where the solution is typically only continuous in $(t\to \mathscr{H})$, with $\mathscr{H}=L^2(\Omega )$, where boundary conditions do not play a role.

REWRITING THE SPATIALLY DISCRETIZED EQUATIONS

In this section, we consider the spatially discretized equations introduced in the section Spatially discretized equations. For the sake of illustration, we assume that we only have Dirichlet boundary conditions. This is not a restriction, since one can always split the boundaries and consider the parts separately.

Direct assignment of the boundary dofs

We now illustrate that the immediate way of assigning the dofs at the boundary, as it is commonly done for inhomogeneous Dirichlet conditions for stationary problems [21], does not simply lead to a system of the form (1).

Consider the linear formulation (6) of Problem the 2.4 with the assignment of the boundary conditions as in Equation (13b):

(14a)$$M_I\dot{v}-A_I\upsilon =f,$$

(14b)$$G\upsilon ={\upsilon }_{\Gamma }=u,$$

(14c)$$\upsilon (0)=\alpha .$$

Then, with the partitioning of M_I and A_I as in Equation (9), the state equation reads:

$$\left[\begin{array}{cc}{M}_{II}& M_{I\Gamma }\end{array}\right]\left[\begin{array}{c}{\dot{\upsilon }}_I\\ {\dot{\upsilon }}_{\Gamma }\end{array}\right]=A_{II}{\upsilon }_I+A_{I\Gamma }{\upsilon }_{\Gamma }+f$$

which, having inserted Equation (14b), gives:

(15)$$M_{II}{\dot{v}}_I=A_{II}{v}_I+f+A_{I\Gamma }u-M_{I\Gamma }\dot{u}.$$

System (15) is not of the form (2) because of the appearance of $\dot{u}$.

Remark 3.1. One can define a new input as $\tilde{u}:=\dot{u}$ and consider the system:

$$\left[\begin{array}{cc}1& 0\\ 0& M_{II}\end{array}\right]\frac{d}{dt}\left[\begin{array}{c}u\\ {\upsilon }_I\end{array}\right]-\left[\begin{array}{c}0\\ A_{II}({\upsilon }_I+u)+f\end{array}\right]=\left[\begin{array}{c}1\\ -M_{I\Gamma }\end{array}\right]\tilde{u}.$$

This approach uses a so-called dynamical controller that is defined via a differential relation. As pointed out in Ref. [22], for a dynamical controller one can set the initial value to zero to circumvent the expected inconsistencies mentioned in Remark 2.6.

Although Rothe’s method is out of consideration, since it leads to a sequence of algebraic equations rather than to a differential equation, it is worth mentioning that it implicitly approximates the time-derivative of the Dirichlet conditions as they appear in Equation (15).

Remark 3.2. Using Rothe’s method of discretizing Problem (2) in time, first, e.g., via using an explicit Euler scheme with a time step length $\tau$, and subsequently discretizing in space, leads to systems of type:

$$M_I{\upsilon }^{k+1}=M_I{\upsilon }^k+\tau (A_I{\upsilon }^k+f^k),$$

$$G{\upsilon }^{k+1}={\upsilon }_{\Gamma }^{k+1}=u^{k+1},$$

where the superscript $k$ relates to the values of the functions at the kth time instance. If at the previous time instance ${\upsilon }_{\Gamma }^k=u^k$, then the direct assignment at the current instance gives:

$$M_{II}{\upsilon }_I^{k+1}=M_{II}{\upsilon }_I^k+\tau (A_{II}{\upsilon }_I^k+f^k)+\tau A_{I\Gamma }{u}^k-M_{I\Gamma }[u^{k+1}-u^k],$$

which can be seen as time-discrete approximation to Equation (15).

Lifting of the boundary conditions

These approaches base on a lifting $\tilde{\upsilon }$ that fulfills the boundary conditions for all time and the decoupling of the solution $\upsilon =y+\tilde{\upsilon }$, see [23] for an example with linearized Navier–Stokes equations.

We consider the linear time-invariant case (7) and assume that f = 0.

At time $t\in [0,T]$, we define a lifting as:

(16)$$\tilde{\upsilon }(t)=\left[\begin{array}{c}{\tilde{\upsilon }}_I(t)\\ u(t)\end{array}\right].$$

Then, considering Equation (8) with $\upsilon =y+\tilde{\upsilon }$ and splitting A_I and M_I as in Equation (9), we find that $y_{\Gamma }=0$ and we obtain the relation:

$$M_{II}{\dot{y}}_I=A_{II}{y}_I+A_I\tilde{\upsilon }-M_I\dot{\tilde{\upsilon }},y_I(0)={\alpha }_I-{\tilde{\upsilon }}_I(0).$$

We use the abbreviation ${\overline{A}}_{II}=M_{II}^{-1}{A}_{II}$ and, with the well-known solution representation, we obtain that:

$$y_I(t)=e^{{\overline{A}}_{II}t}({\alpha }_I-{\tilde{\upsilon }}_I(0))+{\displaystyle {\int }_0^t{e}^{{\overline{A}}_{II}(t-s)}}{M}_{II}^{-1}(A_I\tilde{\upsilon }-M_I\dot{\tilde{\upsilon }}(s))\text{d}s.$$

After an integration by parts, we find that:

$$y_I(t)=e^{{\overline{A}}_{II}t}({\alpha }_I-{\tilde{\upsilon }}_I(0))+{\displaystyle {\int }_0^t{e}^{{\overline{A}}_{II}(t-s)}}(A_I\tilde{\upsilon }(s)-{\overline{A}}_{II}{M}_{II}^{-1}{M}_I\tilde{\upsilon }(s))\text{d}s$$

$$-M_{II}^{-1}{M}_I\tilde{\upsilon }(t)+e^{{\overline{A}}_{II}t}{M}_{II}^{-1}M\tilde{\upsilon }(0).$$

Using that $M_I\tilde{\upsilon }=M_{II}{\tilde{\upsilon }}_I+M_{I\Gamma }u$ and having regrouped the terms, we conclude that ${\widehat{\upsilon }}_I:=y_I+M_{II}^{-1}M\tilde{\upsilon }={\upsilon }_I+M_{II}^{-1}{M}_{I\Gamma }u$ fulfills the ordinary differential equation (ODE):

(17)$$M_{II}{\dot{\widehat{\upsilon }}}_I=A_{II}{\widehat{\upsilon }}_I+Bu,{\widehat{\upsilon }}_I(0)={\alpha }_I+M_{II}^{-1}{M}_{I\Gamma }u(0),$$

with

$$B=[A_{II}{M}_I^{-1}{M}_{I\Gamma }-A_{I\Gamma }].$$

The actual solution is easily retrieved from $\widehat{\upsilon }={\upsilon }_I+M_{II}^{-1}{M}_{I\Gamma }u$. Note, however, the dependency of the initial value on $u(0)$ in Equation (17).

Remark 3.3. The dependency of the initial value on $u(0)$ is due to the ansatz that assumes smoothness of $\tilde{\upsilon }$, which then extends to the boundary nodes. Accordingly, at the boundary, the initial value needs to be consistent with the control at time t = 0, cf. Remark 2.6.

This is not an issue in practical applications where the determination of a control law does not depend on the initial value for the state like it is the case in linear-quadratic optimal control.

Remark 3.4. We find it worth pointing out, that the system (17) does not depend on the choice of the lifting (16) and, thus, includes in particular the lifting by means of the harmonic extension of the boundary values into the inner.

Split mass matrix lifting

For the particular choice of the lifting:

$$\tilde{\upsilon }(t)=\left[\begin{array}{c}-M_{II}^{-1}{M}_{I\Gamma }u(t)\\ u(t)\end{array}\right]$$

which leads to $M_I\dot{\tilde{\upsilon }}(t)=0$ for all time t, the application for nonlinear systems is straightforward. Considering again, $y=\upsilon -\tilde{\upsilon }$, and the nonlinear case of Problem 2.4, one arrives at the ODE:

$$M_{II}{\dot{y}}_I=A_I(y_I+\tilde{\upsilon }(u))+f,y_I(0)={\alpha }_I+M_{II}^{-1}{M}_{I\Gamma }u(0).$$

Again, the actual solution is easily obtained by a backwards substitution ${\upsilon }_I=y_I+{\tilde{\upsilon }}_I=y_I-M_{II}^{-1}{M}_{I\Gamma }u$, but the initial value depends on the possibly unknown input u.

Remark 3.5. A lifting as defined in this chapter leads to an ODE of the desired form. In a forthcoming work, we will investigate similar manipulations on the abstract equations. If the proposed algebraic splitting has a counterpart in infinite dimensions, then one can expect well posedness of the transformed system also for every finer discretizations.

Remark 3.6. For linear time-dependent cases, similar formulas can be derived using fundamental solution matrices or transition matrices. Also, the split mass matrix approach is readily applicable and gives a system of type (2).

Incorporation of the boundary data via Lagrange multiplier

We consider the formulation of Problem 2.5:

(18a)$$M\dot{\upsilon }(t)-A(t,\upsilon (t))-G^{\text{T}}\lambda (t)=f(t),$$

(18b)$$Gv(t)=u(t),$$

with the Lagrangian multiplier $\lambda$.

The saddle point structure is similar to the velocity-pressure formulation of Navier–Stokes equations, where the pressure can be interpreted as the multiplier that couples the divergence constraint to the momentum equation. In particular, it is a special case of semi-explicit index-2 DAEs as they were considered, e.g., in Ref. [24]. Thus, the formulations for the treatment of the boundary conditions that we propose in this section are adaptions of algorithms for the numerical time integration of Navier–Stokes equations or, more general, DAEs of index 2.

Decoupling by projection. In the considered case, $G$ has the form $G=[\begin{array}{cc}0& I\end{array}]$ and M is symmetric positive definite. Thus, we can define:

$$P:=I-M^{-1}{G}^{\text{T}}{S}^{-1}G,S:=G{M}^{-1}{G}^{\text{T}},\text{and}{Q}^{-}:=S^{-1}G{M}^{-1}.$$

With this, system (8) can be equivalently [15] reformulated as $(\upsilon ,\lambda )=({\upsilon }_i+{\upsilon }_g,\lambda )$, where the transformed variables are the solutions of:

(19)$${\upsilon }_g=M^{-1}{G}^{\text{T}}{S}^{-1}u,$$

(20)$$\lambda =-Q^{-}A(v_g+v_i)-Q^{-}f-Q^{-}M{\dot{\upsilon }}_{\Gamma },$$

and

(21)$${\dot{\upsilon }}_i-P{M}^{-1}A({\upsilon }_i+M^{-1}{G}^T{S}^{-1}u)=P{M}^{-1}f.$$

Note that the differential Equation (21) is of type (1).

With MP = P^TM, in the linear case, we can write the differential equation for ${\upsilon }_i$ as:

$$M{\dot{\upsilon }}_i-P^{\text{T}}A{\upsilon }_i=P^{\text{T}}f+P^{\text{T}}Bu,$$

with $B:=A{M}^{-1}{G}^{\text{T}}{S}^{-1}$. In the nonlinear case, the input appears inside the nonlinearity.

Remark 3.7. Since $n_d\ll n_{\upsilon }$, i.e., the number of dofs associated with the boundary is small if compared to the number of inner nodes, an explicit realization of the projection P is feasible. This is different from the analogue for the Navier–Stokes equation, where the dimension of the subspace of the divergence free functions equals the dimension of the pressure space and, thus, can be large.

Remark 3.8. The variable ${\upsilon }_i$ has zero values at the boundary at all time. Thus, if one only considers the ODE (21) for $v_i$, there is no problem of possibly inconsistent initial values due to the chosen control, cf. Remark 2.6. However, a given initial value has to fulfill also (19).

Regularization via penalization. If one adds the term $\alpha \lambda (t)$, $0<\alpha \ll 1$, to the left-hand side of Equation (18b), one can solve for the multiplier and eliminate it from the differential part:

$$M\dot{\upsilon }(t)-A(t,\upsilon (t))+\frac{1}{\alpha }{G}^{\text{T}}G\upsilon =f(t)+\frac{1}{\alpha }{G}^{\text{T}}u.$$

This approach is known as penalty scheme and pressure penalization in the numerical integration of multibody and Navier–Stokes systems, respectively, cf., e.g., [25,26]. The method is straightforward to implement but comes with the need of a proper choice of the penalization parameter. The main difficulty is that a small parameter α not only increases the quality of the approximation of the constraints but also increases the stiffness of the resulting ODE.

INCORPORATION VIA VARIATIONAL FORMULATIONS AND THEIR DISCRETIZATIONS

In its most general form, the variational or weak incorporation of the Dirichlet boundary conditions is derived from Problem 2.1 as follows. Instead of considering the constraint (4b) one adds a penalty term to the variational formulation of the dynamic Equation (4a):

(22)$$\dot{\upsilon }(t)-\mathcal{A}(t,\upsilon (t))+\frac{1}{\alpha }\lambda \prime (\gamma \upsilon (t)-\mathcal{U}(t))=ℱ(t),$$

where $\lambda \prime :\mathcal{Q}\to \mathcal{V}\prime$ and α is a small penalization parameter. Then, for various choices of $\lambda$ and $\mathcal{Q}$, various weak incorporations of the Dirichlet conditions are realized. For example, defining $\lambda \prime$ through:

$${\langle \lambda \prime q,\varphi \rangle }_{\mathcal{V}\prime ,\mathcal{V}}={\displaystyle {\int }_{\Gamma }\text{(}q,\varphi \text{)}}\text{d}\gamma$$

for a $q\in \mathcal{Q}$ and for any $\varphi \in \mathcal{V}$, one obtains the penalized Robin approximation described in the section Penalized Robin in this paper.

Ultra weak formulations

The basic concept of the ultra weak variational formulation is the transfer of smoothness requirements from the test space to the trial space. In numerical experiments, in a conforming discretization, this concept will require special finite element spaces that are not part of common finite element libraries. We will introduce the formulation and a nonconforming discretization suitable for a direct implementation.

Let $\Phi =W^{2,2}(\Omega )\cap W_0^{1,2}(\Omega )$ and consider the diffusion Equation (5) with β = 0. We call $\upsilon$ an ultra weak solution if:

(23)$${\displaystyle {\int }_{\Omega }\text{(}\dot{\upsilon },\varphi \text{)}}\text{d}\omega -{\displaystyle {\int }_{\Omega }\nu }\text{(}\upsilon ,\Delta \varphi ) \text{d}\omega ={\langle f,\varphi \rangle }_{\Phi \prime ,\Phi }-{\displaystyle {\int }_{\Gamma }\nu }\left(g,\frac{\partial \varphi }{\partial n}\right)\text{d}\gamma$$

for all $\varphi \in \Phi$, cf. [8]. The abstract Equation (23) indicates that a spatial discretization may lead to a system in the form of (2). The difficulty with a conforming approach, however, lies in the definition of matching test functions of high regularity with zero boundary values and suitable ansatz functions.

A possible approach is to drop the requirement that the test functions have zero boundary conditions and to consider:

$${\displaystyle {\int }_{\Omega }\text{(}\dot{v},\varphi \text{)}}\text{d}\omega -{\displaystyle {\int }_{\Omega }\nu }\text{(}v,\Delta \varphi )\text{d}\omega -{\displaystyle {\int }_{\Gamma }\nu }\left(\frac{\partial v}{\partial n},\varphi \right)\text{d}\gamma ={\langle f,\varphi \rangle }_{\Phi \prime ,\Phi }-{\displaystyle {\int }_{\Gamma }\nu }\left(g,\frac{\partial \varphi }{\partial n}\right)\text{d}\gamma ,$$

which can be numerically approximated with standard finite element spaces of sufficiently high regularity.

With the nonconforming ansatz spaces $V\subset W_0^{1,2}(\Omega )$, an approximation to (23) that is readily realizable is given through $v\in V$ that satisfies:

(24)$${\displaystyle {\int }_{\Omega }(\dot{\upsilon },\varphi )}\text{d}\omega +{\displaystyle {\int }_{\Omega }\nu }\text{(}\nabla \upsilon ,\nabla \varphi \text{)d}\omega ={\displaystyle {\int }_{\Omega }\text{(}f,\varphi \text{)}}\text{d}\omega -\nu {\displaystyle {\int }_{\Gamma }\left(g,\frac{\partial \varphi }{\partial n}\right)}\text{d}\gamma$$

for all $\varphi \in V$, cf. [27, Ch. 5.2.1]. This approximation of the solution of the boundary value problem through functions with zero boundaries necessarily leads to a solution of L² regularity regardless of possibly higher regularity of the problem. We have observed this low regularity in experiments using explicit schemes for time integration. However, in the reported numerical tests that use implicit schemes, the discretization (24) leads to decent approximations.

The numerical approximation to ultra weak solutions of elliptic problems with boundary control as proposed in Ref. [7] uses a finite dimensional ansatz space $V\subset H^1(\Omega )$ and as the test space $W:=V\cap H_0^1(\Omega )$, see Refs. [8,10] for the extension to parabolic problems. The elliptic case then reads, find $\upsilon \in V$, such that:

(25a)$${\displaystyle {\int }_{\Omega }\nu }\text{(}\nabla v,\nabla \varphi \text{) d}\omega ={\displaystyle {\int }_{\Omega }(f,\varphi )}\text{ d}\omega ,\text{for all }\varphi \in V\cap H_0^1(\Omega ),$$

(25b)$$\upsilon ={\Pi }_V(u)\text{ on }\Gamma ,$$

where ${\Pi }_V$ is the $L^2$ projection in $L^2(\Gamma )$ onto the grid induced by the triangulation that defines V. Using the spaces and formulations of (9) for a spatial discretization of a parabolic problem, one obtains a system that is the same as (7) apart from the appearance of the projector ${\Pi }_V$ in the boundary term (14b). Anyways, the elimination of the boundary nodes will lead to a system like (15) that includes the time derivative $\dot{u}$ of the control. A discontinuous Galerkin ansatz for the time discretization as used in Ref. [10] includes $\dot{u}$ implicitly in the same way as Rothe’s method, cf. Remark 3.2.

Remark 4.1. Since the known numerical approaches that base on (9) do not lead to systems of type (2) or (1), we did not consider them in the numerical experiments in this manuscript. However, the lifting, cf. the section Lifting of the boundary conditions, and the projection approach, cf. the section Decoupling by projection, readily apply to the formulation of the boundary term that includes the projector ${\Pi }_V$. The inclusion of the projection is necessary for well posedness of the Dirichlet control problem for the case of less regular boundary controls [7,8,10].

Nitsche variational formulation

A variant of the standard weak formulation of the pure diffusion, cf. (2) with β = 0, as proposed in Ref. [28] for the stationary Poisson equation reads:

(26)$${\displaystyle {\int }_{\Omega }(\dot{\upsilon },\varphi )}\text{d}\omega +{\displaystyle {\int }_{\Omega }\nu }\text{(}\nabla \upsilon ,\nabla \varphi ) \text{d}\omega -{\displaystyle {\int }_{\Gamma }\nu }\left(\frac{\partial v}{\partial n},\varphi \right)\text{d}\gamma -{\displaystyle {\int }_{\Gamma }\nu }\left(\upsilon ,\frac{\partial \varphi }{\partial n}\right) \text{d}\gamma +c_{\gamma }{\displaystyle {\int }_{\Gamma }(\upsilon ,\varphi )}\text{d}\gamma ={\langle f,\varphi \rangle }_{\Phi \prime ,\Phi }-{\displaystyle {\int }_{\Gamma }\nu }\left(g,\frac{\partial \varphi }{\partial n}\right)\text{d}\gamma +c_{\gamma }{\displaystyle {\int }_{\Gamma }(g,\varphi )}\text{d}\gamma$$

for all $\varphi \in \Phi =W^{1,2}(\Omega )$. The formulation is derived by considering the cost functional:

$$\mathcal{J}(w)={\displaystyle {\int }_{\Omega }\nu }\text{(}\nabla w,\nabla w)\text{d}\omega -2{\displaystyle {\int }_{\Gamma }\nu }\left(\frac{\partial w}{\partial n},w\right)\text{d}\gamma +c_{\gamma }{\displaystyle {\int }_{\Gamma }(w,w)}\text{ d}\gamma ,$$

with a parameter $c_{\gamma }$ and the first order optimality conditions for $\mathcal{J}(w-v)\to \min$, where $\upsilon$ is the solution of the stationary Poisson problem with nonhomogeneous Dirichlet boundary conditions. If for a given mesh $c_{\gamma }$ is chosen sufficiently large, namely $c_{\gamma }\approx h^{-1}$ where $h$ is a characteristic length of the triangulation, then the discretized optimization problem is convex [28, Equation (12)].

For (26), a standard discrete formulation leads to an equation of type (2) with A and B explicitly given, see Ref. [29]. Cf. also [27, Ch. 5.2.2] where nonzero boundary values of $y$ have been assumed.

Penalized Robin

If one approximates the Dirichlet conditions by a Robin-type condition:

$$\upsilon \approx \alpha \frac{\partial \upsilon }{\partial n}+\upsilon =g\text{or}\frac{\partial \upsilon }{\partial n}\approx \frac{1}{\alpha }(g-\upsilon )\text{ on }\Gamma ,$$

with a parameter α that is intended to go to zero, then the boundary conditions are incorporated naturally in the weak formulation of the convection–diffusion operator (6) and a standard finite element discretization leads to a system of type (1). For the pure diffusion case, convergence of the solutions to the actual solution for $\alpha \to 0$ has been shown in several contexts, cf. [12] and the references therein.

NUMERICAL TESTS

We consider two-dimensional convection–diffusion-reaction problems. All setups are directed to actuation at the boundary. In particular, there are no source terms. This excludes the method of manufactured solutions for consistency and convergence checks, where one constructs a solution and derives the corresponding source term and boundary data. In any case, the method of manufactured solution seems not well suited to test the modeling of boundary actuation, since the numerical solution will depend almost exclusively on the volume force; see the test case at the end of this section.

Hence, in order to evaluate the convergence numerically, we compute a reference solution using the naive approach (15) of directly assigning the boundary nodes and a very fine grid in space and time.

We refer to the tested schemes as follows:

• dias – direct assignment of the boundary values – cf. the section Direct assignment of the boundary dofs

• lift – lifting of the boundary conditions via split mass matrix – cf. the section Lifting of the boundary conditions

• proj – incorporation of the constraint via Lagrange multiplier and projections – cf. the section Decoupling by projection

• pena – penalization of the constraint – cf. the section Regularization via penalization

• ncul – nonconforming approximation of ultra week solutions – cf. the section Ultra weak formulations

• nits – approximation via the Nitsche variational formulation – cf. the section Nitsche variational formulation

• pero – relaxation via Robin approximation – cf. the section Penalized Robin

For all test setups, we will check the convergence of dias and that the theoretically equivalent formulations lift and proj give the same results. Also, we will investigate how the relaxed methods pena, nits, and pero perform for different choices of the penalization parameter and for inexact solves of the resulting linear systems. Furthermore, we will investigate how the schemes perform if an iterative solver is applied.

Test setups

We consider several convection–diffusion setups on a two-dimensional square domain. Let $\Omega =[-1,1]\times [-1,1]\subset {ℝ}^2$ be the computational domain with the spatial coordinates $x=(x_0,x_1)$. Let Γ be the boundary with parts Γ₀ to Γ₃ as depicted in Figure 1. All setups model the evolution in time and space of a scalar quantity ρ due to a convection wind β and diffusion with a diffusion coefficient $\nu$, cf. Problem 2.2.

Figure 1.

Illustration of the domain, the arrangement of the boundary segments, a triangulation with N_h = 12, and a snapshot of an approximation to the internal convection–diffusion as described in Test Case 1 at time t = 3.0.

The quantity ρ is seeded into the domain at Γ₀, where we enforce the time-dependent Dirichlet conditions:

(27)$$\rho {|}_{\Gamma }=g(x)u(t):=\frac{1}{2}\left(\sin \left(\pi x_0+\frac{\pi }{2}\right)+1\right)(\cos (2t+\pi )+1).$$

Here, $g(x):=\frac{1}{2}\left(\sin \left(\pi x_0+\frac{\pi }{2}\right)+1\right)$ is the spatial shape function and $u(t):=\cos (2t+\pi )+1$ is the scalar control function. At the remainder boundaries, Γ₁, Γ₂, and Γ₃, depending on the setup, homogeneous Dirichlet or homogeneous Neumann boundary conditions are applied. As the initial value, we set ρ(0) = 0, which is consistent with the control action at time t = 0.

As the first test case, we consider a setup with no convection at the boundary, so that the boundary control is propagated into the domain only by diffusion.

Test case 1 (internal convection–diffusion). Given a convection wind and a diffusion parameter:

$${\beta }_0(x)=\left[\begin{array}{c}-x_1{(x_0-1)}^2{(x_0+1)}^2(x_1-1)(x_1+1)\\ x_0(x_0-1)(x_0+1){(x_1-1)}^2{(x_2+1)}^2\end{array}\right]\text{and}{\nu }_0=0.1,$$

find approximations to the scalar function ρ satisfying:

(28a)$$\dot{\rho }(t)+{\beta }_0·\nabla \rho (t)-{\nu }_0\Delta \rho (t)=0,$$

(28b)$$\rho {|}_{{\Gamma }_0}(t)=gu(t),$$

(28c)$$\rho {|}_{{\Gamma }_1\cup {\Gamma }_2{\Gamma }_3}(t)=0,$$

(28d)$$\rho (0)=0,$$

on given discretizations of the spatial domain $\Omega ={[-1,1]}^2$ and of the time interval $[0,4]$.

As a second test case, we consider a convection–diffusion problem with inflow and outflow, for which the boundary values are also transported into the domain via convection. See Figure 2a for an illustration of the setup.

Figure 2.

Illustration of the distribution of the scalar ρ seeded at the upper boundary after diffusion and convection (a) and additional reaction (b) as described in Test Cases 2 and 3 for N_h = 15 at time t = 3.0.

Test Case 2 (convection–diffusion). Given a convection wind and a diffusion parameter

$${\beta }_1(x)=\frac{1}{10}\left[\begin{array}{c}{x}_0+1\\ -(x_1+1)\end{array}\right]\text{and}{\nu }_1=0.1,$$

find approximations to the scalar function ρ satisfying:

$$\dot{\rho }(t)+{\beta }_1·\nabla \rho (t)-{\nu }_1\Delta \rho (t)=0,$$

$$\rho {|}_{{\Gamma }_0}(t)=gu(t),$$

$$\rho {|}_{{\Gamma }_1\cup {\Gamma }_2}(t)=0,$$

$$\frac{\partial \rho }{\partial \nu }{|}_{{\Gamma }_3}(t)=0,$$

$$\rho (0)=0,$$

on given discretizations of the spatial domain $\Omega ={[-1,1]}^2$ and of the time interval $[0,0.2]$.

The third test case is the same as the second but with an additional reaction source term $r(\rho )=\rho (1-\rho )$ in the dynamical equation. This source term $r$ is positive for values of $0\le \rho \le 1$ and negative elsewhere. Thus, for values of ρ > 0 the reaction pushes ρ towards ρ = 1, cf. Figure 2b.

The considered system, for $t\in (0,1]$, now reads

Test Case 3. Given the wind and the diffusion parameter defined in Test Case 2, find approximations to the scalar function ρ satisfying:

$$\dot{\rho }(t)+{\beta }_1·\nabla \rho (t)-{\nu }_1\Delta \rho (t)=\rho (t)(1-\rho (t)),$$

$$\rho {|}_{{\Gamma }_0}(t)=gu(t),$$

$$\rho {|}_{{\Gamma }_1\cup {\Gamma }_2}(t)=0,$$

$$\frac{\partial \rho }{\partial \nu }{|}_{{\Gamma }_3}(t)=0,$$

$$\rho (0)=0,$$

on given discretizations of the spatial domain $\Omega ={[-1,1]}^2$ and of the time interval $[0,0.2]$.

For all test cases, the spatial discretization is done on a uniform criss-cross triangulation described by the parameter $N_h=\frac{2}{h}$ which is the length of the boundary parts divided by the length of the longest edge of the triangles, see Figure 1. For the discrete function space, we use the parameter $cg$, denoting the polynomial degree of the chosen standard Lagrange elements. For the time discretization, we use a uniform grid of size $N_{\tau }\approx \frac{1}{\tau }$ corresponding to the ratio of the length of the time interval versus the length of one time step. Here and in the following examples, already for the coarsest discretization, the local Peclet number $Pe:=\Vert \beta \left(t\right)\Vert h/v$ is smaller than 1. Thus, we can expect reliable approximations without additional, e.g., upwind, stabilization [30].

For the spatial discretization, we use the Python interface dolfin [31] to the finite element software suite Fenics [32]. Our investigation focusses on the space discretization errors but we will make sure that the time integration error is sufficiently small. For the linear cases, the time integration is done by means of the implicit trapezoidal rule. The nonlinear case is treated implicitly in the linear part and with the Method of Heun in the nonlinear part. The norms are computed using the piecewise trapezoidal rule for the time integration and dolfin’s built-in function error norm that evaluates the $L^2$ norm in the discrete function spaces. In general, we solve the occurring linear equation systems via a direct solver that makes use of the python module scipy’s built-in sparse LU decomposition method. In some tests, we employ the generalized minimal residual method (GMRES) method using the implementation of the python module krypy [33]. The code used can be obtained from the author’s public git repository [34].

By ${\rho }_{h{N}_h,\tau N_{\tau }}^{pcg}$, we denote the approximation to the solution of (10) with the discretization parameters $N_h$, $N_{\tau }$, and $cg$. By $e_{h{N}_h,\tau N_{\tau }}^{pcg}$, we denote the approximation error

$$e_{h{N}_h,\tau N_{\tau }}^{pcg}:={\rho }_{h{N}_h,\tau N_{\tau }}^{pcg}-{\rho }_{\text{ref}}$$

measured in a numerical approximation of the $L^2(0,1;L^2({[-1,1]}^2))$ norm, where ρ_ref is a reference computed with the $cg=2$ scheme with $N_{\tau }=240$ and $N_h=96$.

Convergence tests

In Tables 1 and 2, we list the approximation errors for dias for increasingly fine space and time discretizations for Test Cases 1 and 2. One can see, that the spatial discretization error is dominating, i.e., convergence in the time discretization is only observed for larger values of $N_{\tau }$. This justifies the choice of $N_{\tau }:=240$ as the reference discretization for further error comparisons.

Table 1.

(Time space convergence of dias for linear elements, cf. the section Convergence tests) The approximation error $e_{hNh,\tau N\tau }^{pcg}$ with ρ_ref = ${\text{ρ}}_{h96,\tau 120}^{p2}$ scaled by the inverse of $e_{h6,\tau 30}^{p1}=1.119·{10}^{-1}$ (top) and $e_{h6,\tau 60}^{p2}=3.201·{10}^{-2}$ (bottom) for varying space and time discretizations and for polynomial degree $cg=1$ (top) and $cg=2$ (bottom) for Test Case 1. Cf. also Figure 3a and b illustrating the convergence in space for the finest time discretization (the rightmost columns in the tables).

${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	30	60	120
6	1.0000	1.0026	1.0033
12	0.2608	0.2641	0.2651
24	0.0654	0.0661	0.0671
48	0.0244	0.0163	0.0166
96	0.0215	0.0059	0.0041
${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	60	120	240
6	1.0000	0.9982	0.9978
12	0.2295	0.2272	0.2269
24	0.0482	0.0424	0.0419
48	0.0201	0.0077	0.0063

Table 2.

(Time space convergence of dias for quadratic elements, cf. the section Convergence tests) the approximation error $e_{hNh,\tau N\tau }^{pcg}$ with ρ_ref = ${\text{ρ}}_{h96,\tau 120}^{p2}$ scaled by the inverse of $e_{h6,\tau 30}^{p1}=4.349·{10}^{-4}$ (top) and $e_{h6,\tau 60}^{p2}=8.551·{10}^{-05}$ (bottom) for varying space and time discretizations and for polynomial degree $cg=1$ (top) and $cg=2$ (bottom) for Test Case 2. Cf. also Figure 3c and d illustrating the convergence in space for the finest time discretization (the rightmost columns in the tables).

${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	30	60	120
6	1.0000	0.9997	0.9996
12	0.3696	0.3695	0.3694
24	0.1060	0.1059	0.1059
48	0.0276	0.0275	0.0275
96	0.0071	0.0070	0.0069
${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	60	120	240
6	1.0000	1.0000	0.9999
12	0.1699	0.1699	0.1699
24	0.0316	0.0330	0.0305
48	0.0085	0.0071	0.0071

The errors $e_{hNh,\tau 120}^{pcg}$ for a fixed time discretization and varying space discretizations are plotted in Figure 3 for all three test cases. From the plots, one can see that the equivalent formulations lift and proj give the same results and one can read off the numerically estimated order of spatial convergence (EOC). For the linear elements ($cg=1$), one obtains $EOC=2$ and for the quadratic elements ($cg=2$), one obtains $EOC=2.5$ at a lower error level. The observed order of convergence is not optimal as laid out in the section Convergence tests with volume forcing. For the linear elements, also ncul converges quadratically although with an error that is slightly larger than the one reported for proj and lift.

Figure 3.

Convergence tests for the consistent implementations, cf. the section Convergence tests. The error $e_{h{N}_h,\tau 120}^{pcg}$ for varying space discretizations N_h and for linear (left) and quadratic (right) shape functions. The first row of plots (a and b) corresponds to Test Case 1, the middle row (c and d) to Test Case 2, and the bottom line (e and f) to Test Case 3. The dashed lines indicate the slope of a quadratic convergence the dotted lines indicate a convergence of order 2.5.

For piecewise quadratic shape functions, the scheme ncul delivered good approximations but its convergence rate was estimated as $EOC=0.5$, see Figure 3 (right column), which is much less than expected from theory. A possible explanation for this breakdown is the oscillation that occurs when approximating a step function by a quadratic polynomial. Recall that the scheme ncul enforces a zero value at the boundary, while in the inner it approximates a solution which is not zero at the boundary. The inevitable jump in the solution approximation is seemingly well captured by linear but not by quadratic elements.

Parameter studies for the penalty schemes

For the schemes pena, nits, and pero that depend on a parameter, we investigate the accuracy of the approximation versus the choice of the penalization parameter α, where we have defined the relation $c_{\gamma }=\frac{\nu }{\alpha }$ to fit in Nitsche’s method (26). Judging from the results depicted in Figure 4, for large penalization parameters, the approximation is bad, while for small parameters the accuracy of the consistent approximations is obtained. The Nitsche method nits did not lead to reasonable approximations for large values of a.

Figure 4.

Penalization parameter study, cf. the section Parameter studies for the penalty schemes. The error $e_{h{N}_h,\tau 120}^{p1}$ for varying space discretizations N_h versus the penalization parameter α for the schemes pena (left), pero (middle), and nits (right). The first row of plots (a–c) corresponds to Test Case 1, the middle row (d–f) to Test Case 2, and the bottom line (g and h) to Test Case 3.

The necessity to properly choose the penalization parameter is evident in the errors that are reported for inexact solutions of the resulting linear systems. If one applies GMRES preconditioned with the inverse of the mass matrix, to solve the algebraic equations in every time step, the approximation error of the penalization schemes increases with smaller penalization parameters α. The plots in Figure 5 show this phenomenon. For this investigation, we allowed a relative residual of at most $tol={10}^{-5}$, which is already less than the overall error which is of magnitude 10⁻⁴.

Figure 5.

Penalty schemes and inexact solves, cf. the section Parameter studies for the penalty schemes. The error $e_{h{N}_h,\tau 120}^{p1,tol1e-5}$ for varying space discretizations $N_h$ versus the penalization parameter α for the schemes $pena$ (left), $pero$ (middle), and $nits$ (right), where the occurring algebraic equations are solved via GMRES up to an residual of ${10}^{-5}$. The first row of plots (a–c) corresponds to Test Case 1, the middle row (d–f) to Test Case 2, and the bottom line (g–i) to Test Case 3.

The increase in the approximation error is mainly due to the increase of the magnitude of the right-hand side that scales with $\frac{1}{\alpha }$. In fact, having solved the exemplary linear system $Ax=f$ up to a relative residual of $tol$, one has that:

$$\frac{\Vert Ax-f\Vert }{\Vert f\Vert }=\text{tol}\text{or}\Vert Ax-f\Vert =\text{tol}\Vert f\Vert$$

which means that for larger right-hand sides f, the absolute residual $\Vert Ax-f\Vert$ can be larger. A remedy is to control the absolute residual which can be done by correcting the provided relative residual by a factor $\text{tolcor}\text{=}\text{min}\left\{\frac{\text{1}}{\Vert f\Vert },1\right\}$, where $f$ denotes the current right-hand side. In Figure 6d–f, we have reported the discrete $L^2(0,0.2)$ norm of $tolcor$ for Test Case 2. Applying this correction, that scales with $\frac{1}{\alpha }$, one recovers the approximation properties of exact solves over the whole range of α, cf. Figures 6g– i and 4d–f.

Figure 6.

Penalty schemes and absolute tolerances, cf. the section Parameter studies for the penalty schemes. The error $e_{h{N}_h,\tau 120}^{p2,tol1e-5}$ for a fixed relative residual $tol={10}^{-5}$ (top row), the correction of the residual tolcor (middle row), and for a fixed absolute residual $abstol={10}^{-5}$ (bottom row) for varying space discretizations $N_h$ versus the penalization parameter α for the schemes pena (left), pero (middle), and nits (left) for Test Case 2.

GMRES performance

In this test setup, we investigated how the different but mainly equivalent formulations of the same problem affect the performance of an iterative solver. Therefore, we fixed the time and space discretization $N_h=48$ and $N_{\tau }=120$ and, for polynomial degrees $\mathrm{cg}=1,2$, we considered the simulations of Test Case 1, 2, and 3 if the resulting linear equations are solved using GMRES up to a relative residual smaller than $tol={10}^{-7}$. The results are listed in Tables 3 and 4.

Table 3.

Performance of GMRES within various formulations, see section GMRES performance. The averaged number of iterations per time-step $av.\#its$ and the approximation error for several methods and all three test cases for $N_h=48$, $N_{\tau }=120$, and linear elements ($cg=1$ ) in the case that the resulting linear equations are solved using GMRES up to a relative residual of $tol={10}^{-7}$.

	Test Case 1		Test Case 2		Test Case 3
	$\mathit{a}\mathit{v}.\mathit{\#}\mathit{i}\mathit{t}\mathit{s}$	${\mathit{e}}_{\mathit{h}48,\mathit{\tau }120}^{\mathit{p}1,\mathit{t}ol1\mathit{e}-7}$	$\mathit{a}\mathit{v}.\mathit{\#}\mathit{i}\mathit{t}\mathit{s}$	${\mathit{e}}_{\mathit{h}48,\mathit{\tau }120}^{\mathit{p}1,\mathit{t}ol1\mathit{e}-7}$	$\mathit{a}\mathit{v}.\mathit{\#}\mathit{i}\mathit{t}\mathit{s}$	${\mathit{e}}_{\mathit{h}48,\mathit{\tau }120}^{\mathit{p}1,\mathit{t}ol1\mathit{e}-7}$
proj	41.7	1.9·10−3	10.6	1.2·10−5	10.6	1.1·10−5
lift	41.7	1.9·10−3	10.6	1.2·10−5	10.6	1.1·10−5
pero	60.8	4.0·10−3	14.2	7.8·10−6	14.2	7.6·10−6
pena	48.9	9.7·10−3	15.5	1.2·10−5	15.5	1.1·10−5
ncul	43.2	3.7·10−3	10.6	1.1·10−5	10.6	1.1·10−5

The colored cells contain the lowest measured values.

Table 4.

(Performance of GMRES within various formulations, see section GMRES performance) The averaged number of iterations per time-step $av.\#its$ and the approximation error for several methods and all three test cases for $N_h=48$, $N_{\tau }=120$, and quadratic elements ($cg=2$ ) in the case that the resulting linear equations are solved using GMRES up to a relative residual of $tol={10}^{-7}$.

	Test Case 1		Test Case 2		Test Case 3
	$\mathit{a}\mathit{v}.\mathit{\#}\mathit{i}\mathit{t}\mathit{s}$	${\mathit{e}}_{\mathit{h}48,\mathit{\tau }120}^{\mathit{p}2,\mathit{t}ol1\mathit{e}-7}$	$\mathit{a}\mathit{v}.\mathit{\#}\mathit{i}\mathit{t}\mathit{s}$	${\mathit{e}}_{\mathit{h}48,\mathit{\tau }120}^{\mathit{p}2,\mathit{t}ol1\mathit{e}-7}$	$\mathit{a}\mathit{v}.\mathit{\#}\mathit{i}\mathit{t}\mathit{s}$	${\mathit{e}}_{\mathit{h}48,\mathit{\tau }120}^{\mathit{p}2,\mathit{t}ol1\mathit{e}-7}$
proj	83.4	2.5·10−4	20.1	6.3·10−7	20.1	6.1·10−7
lift	83.4	2.5·10−4	20.1	6.0·10−7	20.1	6.1·10−7
pero	106.7	4.6·10−3	24.9	1.1·10−5	24.9	1.1·10−5
pena	71.2	3.8·10−2	20.6	4.2·10−6	20.6	4.4·10−6
ncul	84.9	7.7·10−2	20.3	1.1·10−4	20.3	1.1·10−4

The colored cells contain the lowest measured values.

The residuals were calculated in the inner product induced by the inverse of the mass matrices which was achieved by using the inverses as a preconditioner. At each timestep, as initial guesses, we took the values obtained by linear extrapolation on the bases of the two latest values. The parameter α was set to $\alpha =1$ for pena and $\alpha ={10}^{-3}$ for pero which corresponds to the optimal values for the $cg=1$ case, cf. Figure 5.

As a performance measure that is comparatively independent of the sophistication of the implementation, we took the averaged numbers of iteration per timestep that were needed to obtain a residual below $tol={10}^{-7}$. A second quality measure was the resulting approximation error with respect to the reference solution. In all tests, the methods proj and lift took the least number of iterations. In some cases, in terms of approximation quality, they were outperformed by pero, but at the price of significantly more necessary iterations. The scheme ncul performs similar to proj and lift for $cg=1$. For $cg=2$ the approximation was much worse as it was already observed in Ref. [22]. At almost all tests, the penalization schemes needed more iterations and lead to worse approximations if compared to the consistent schemes. Note, however, that the choice of the penalization parameters was certainly not optimal for the $cg=2$ cases.

Convergence tests with volume forcing

In the beginning of the section Numerical tests, we have mentioned that the method of manufactured solutions is not suitable for boundary controlled processes. This is intuitively clear since for every finer discretizations the weight of a boundary tends to zero if compared to a surface or volume patch. More concretely, in two spatial dimensions, the number of nodes at the boundary grows linearly, while the number of nodes in the inner grows at least quadratically. Thus, if the boundary conditions are merely an extension of a volume force, the volume force will dominate over what happens at the boundary.

To back this assertion by a numerical experiment, we consider Test Cases 1 and 2 (see the section Test setups) but with an additional volume force in Equation (10a) corresponding to the constructed solution:

$${\rho }_{\text{ref}}=\frac{1}{8}\left(\sin \left(x_0\pi +\frac{\pi }{2}\right)+1\right)\left(\sin \left(\frac{x_1}{2}\pi \right)+1\right)(1+x_1)u(t),$$

with $u$ as in Equation (27). The solution ρ_ref is constructed such that at ${\Gamma }_0$ it coincides with the boundary control function defined in Equation (27) and such that it is zero at the remaining boundaries. Also, it holds that $\frac{\partial {\rho }_{\text{ref}}}{\partial \nu }{|}_{{\Gamma }_3}=0$ as required for the setup of Test Case 2.

Taking the method lift and tabulating the approximation errors for varying time and space discretization, for linear elements, we find spatial convergence orders $EOC=2$, i.e., doubling $N_h$ reduces the error by a factor of $2^{-2}$. For quadratic elements, we find $EOC=3$, i.e., doubling $N_h$ reduces the error by a factor of $2^{-3}$, cf. Tables 5 and 6, and Figure 7. The convergence order is as expected for stationary problems and, for the quadratic ansatz functions, significantly better than in the previous experiments, cf., in particular, Table 1 and Figure 3b and d. This indicates that the boundary conditions are not optimally considered by standard discretization schemes. Moreover, this insufficiency is not captured by numerical tests with systems that are driven by a volume force.

Table 5.

(Time space convergence of lift with volume forcing, cf. section Convergence tests) The approximation error $e_{hNh,\tau N\tau }^{p1}$ scaled by the inverse of $e_{h6,\tau 30}^{p1}=9.7149·{10}^{-2}$ for linear ansatz functions (top) and $e_{hNh,\tau N\tau }^{p1}$ scaled by the inverse of $e_{h6,\tau 60}^{p2}=5.288·{10}^{-3}$ for quadratic ansatz functions (bottom) with ρ_ref explicitly given for varying space and time discretizations for Test Case 1.

${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	30	60	120
6	1.0000	0.9975	0.9975
12	0.2720	0.2579	0.2579
24	0.1064	0.0652	0.0651
48	0.0797	0.0172	0.0163
96	0.0766	0.0067	0.0041
${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	60	240	960
6	1.0000	0.9429	0.8810
12	0.3681	0.1049	0.1018
24	0.3488	0.0258	0.0124
48	0.3485	0.0218	0.0021

Table 6.

(Time space convergence of lift with volume forcing, cf. section Convergence tests) The approximation error $e_{hNh,\tau N\tau }^{p1}$ scaled by the inverse of $e_{h6,\tau 30}^{p1}=1.29·{10}^{-4}$ for linear ansatz functions (top) and $e_{hNh,\tau N\tau }^{p1}$ scaled by the inverse of $e_{h6,\tau 60}^{p2}=7.234·{10}^{-6}$ for quadratic ansatz functions (bottom) with ρ_ref explicitly given for varying space and time discretizations for Test Case 2.

${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	30	60	120
6	1.0000	1.0773	0.9992
12	0.2744	0.2740	0.2740
24	0.0614	0.0610	0.0610
48	0.0153	0.0152	0.0152
96	0.0039	0.0038	0.0038
${\mathit{N}}_{\mathit{h}}\backslash {\mathit{N}}_{\mathit{\tau }}$	60	120	240
6	1.0000	0.9998	0.9997
12	0.1175	0.1174	0.1174
24	0.0140	0.0139	0.0139
48	0.0022	0.0017	0.0017

Figure 7.

Spatial convergence with manufactured solutions, cf. section Convergence tests. The error $e_{h{N}_h,\tau N_{\tau }}^{pcg}$ for Test Case 1 (a) for Test Case 2 (b) for sufficiently fine time discretizations $N_{\tau }$, for varying space discretizations $N_h$, and for linear and quadratic shape functions. The dashed lines indicate the slope of a quadratic convergence the dotted lines indicate a convergence of order 3.

CONCLUSION

We have listed common numerical schemes and introduced a projection-based method for problems with time-dependent Dirichlet boundary conditions. We have made the distinction between consistent schemes and relaxed schemes that depend on a penalization parameter.

Using a reference solution on a fine discretization, we investigated the order of convergence of the space discretization for the different schemes. The estimated order of convergence was in between $EOC=2$ and $EOC=2.5$ which is not satisfactory. Similar tests but with a volume force led to an $EOC=3$, the quadratic elements. This result suggests that boundary-driven problems are not treated optimally in the considered finite element schemes. A numerical analysis would be needed to detect the source of the breakdown and to find remedies like, maybe, the boundary concentrated Finite Element approximation [35]. Apart from that, the results as a whole show that the method of manufactured solutions is not well suited for the numerical investigation of spatial convergence of boundary actuation-driven setups.

The relaxed schemes showed the same accuracy as the consistent schemes, but only at certain ranges of the penalization parameter value. If one solves the algebraic equations with high accuracy, one only has to choose the penalization small enough. However, if the algebraic equations are solved iteratively up to a certain residual, then the approximation gets worse again for smaller penalization parameters. This effect might be partially due to an ill-conditioning of the system which might be cured by a suitable preconditioner. The main factor, however, is that for small penalization parameters α the residual is dominated by the penalization term. As a remedy, one can consider absolute residuals as convergence criteria. Conversely, that means that one has to prescribe relative residuals that scale with α which is not practical for small α.

In addition to the approximation quality, we have investigated the performance of GMRES applied within the various schemes. In these tests, as expected, the consistent schemes outperformed the schemes with a penalization.

Based on the results, depending on the situation, we speak out in favor of certain methods as follows. In view of minimal effort for implementation, pero and ncul are the methods of choice. If one wants to invest some time in implementation, proj and lift are better choices since they provide reliable approximations independent of parameters and for higher-order elements and they perform better in iterative schemes. If one can afford the incorporation of the projector, in particular for optimization, proj might be preferable over lift since a possibly inconsistent initial value is not an issue here.

A main motivation of the survey was that standard model reduction or optimal control approaches are readily applicable to systems of distributed type like (2). In a forthcoming paper, we will investigate how well the proposed formulations work in control setups. Also the consistency of the reformulations with the abstract equations is still open and subject to ongoing work.

References

1. King R, editor. Active flow control. Papers contributed to the conference ‘Active flow control 2006’; 2006 Sep 27–29; Berlin (Germany): Springer; 2007.

2. King R, editor. Active flow control II: Papers contributed to the conference ’Active Flow Control II 2010’; 2010 May 26–28. Notes on numerical fluid mechanics and multidisciplinary design. Berlin (Germany): Springer; 2010.

3. , . Boundary control of PDEs. A course on Backstepping designs. Philadelphia (PA): SIAM; 2008.

4. , , , . Representation and control of infinite-dimensional systems. Vol. I. Basel (Switzerland): Birkhäuser; 1992.

5. Tröltzsch F. Optimale Steuerung partieller Differentialgleichungen. Wiesbaden (Germany): Vieweg+Teubner; 2009.

6. Lions JL. Optimal control of systems governed by partial differential equations. Berlin (Germany): Springer; 1971.

7. Berggren M. Approximations of very weak solutions to boundary-value problems. SIAM J Numer Anal. 2004;42(2):860–77.

8. , , , . Optimization with PDE constraints. Dordrecht (Netherlands): Springer; 2009.

9. , , . Error analysis for a finite element approximation of elliptic Dirichlet boundary control problems. SIAM J Control Optim. 2013;51(3):2585–611.

10. , , . A finite element method for Dirichlet boundary control problems governed by parabolic PDEs. ArXiv e-prints, Oct. 2014. Available from: http://arxiv.org/abs/1410.0136

11. , . Constrained Dirichlet boundary control in L² for a class of evolution equations. SIAM J Control Optim. 2007;46(5):1726–53.

12. , , . A penalized Robin approach for solving a parabolic equation with nonsmooth Dirichlet boundary conditions. Asymptotic Anal. 2003;34(2):121–36.

13. , . Error estimates for the numerical approximation of Dirichlet boundary control for semilinear elliptic equations. SIAM J Control Optim. 2006;45(5):1586–611.

14. , . LQG-balanced truncation low-order controller for stabilization of laminar flows. In: King R, editor. Active flow and combustion control 2014, Volume 127 of notes on numerical fluid mechanics and multidisciplinary design. Berlin: Springer; 2015. p. 365–379.

15. Heiland J. Decoupling and optimization of differential-algebraic equations with application in flow control (PhD thesis). Technical University of Berlin; 2014. Available from: http://opus4.kobv.de/opus4-tuberlin/frontdoor/index/index/docId/5243

16. Braess D. Finite Elemente. Theorie, schnelle Löser und Anwendungen in der Elastizitätstheorie. 4th revised and extended edition. Berlin (Germany): Springer; 2007.

17. , . Finite element methods for Navier–Stokes equations theory and algorithms. Berlin (Germany): Springer; 1986.

18. Roubíček T. Nonlinear partial differential equations with applications. Basel (Switzerland): Birkhäuser; 2005.

19. , . Treating inhomogeneous essential boundary conditions in finite element methods and the calculation of boundary stresses. SIAM J Numer Anal. 1992;29(2):390–424.

20. Altmann R. Index reduction for operator differential-algebraic equations in elastodynamics. Z Angew Math Mech. 2013;93(9):648–64.

21. , , , editors. Automated solution of differential equations by the finite element method, Volume 84 of Lect. Notes Comput Sci Eng. 1st edition. Springer-Verlag; 2012.

22. , . Stabilization of parabolic nonlinear systems with finite dimensional feedback or dynamical controllers: application to the Navier–Stokes system. SIAM J Control Optim. 2011;49(2):420–63.

23. Raymond J-P. Stokes and Navier–Stokes equations with nonhomogeneous boundary conditions. Ann Inst H Poincaré Anal Non Linéaire. 2007;24(6):921–51.

24. , , . The numerical solution of differential-algebraic systems by Runge-Kutta methods. Berlin (Germany): Springer; 1989.

25. , . On the stabilizing properties of energy-momentum integrators and coordinate projections for constrained mechanical systems. In: , , , editors. Multibody dynamics, Volume 4 of computational methods in applied sciences. Amsterdam (Netherlands): Springer; 2007. p. 49–67.

26. Rannacher R. On the numerical solution of the incompressible Navier–Stokes equations. Z Angew Math Mech. 1993;73(9):203–16.

27. , . Control theory for partial differential equations: continuous and approximation theories I. Abstract parabolic systems. Cambridge (UK): Cambridge University Press; 2000.

28. Nitsche J. Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind. Abh Math Semin Univ Hambg. 1971;36(1):9–15.

29. Schieweck F. Uniformly stable mixed hp -finite elements on multilevel adaptive grids with hanging nodes. ESAIM Math Model Numer Anal. 2008;42(3):493–505.

30. , . Finite element methods for computational fluid dynamics: a practical guide, Volume 14. Philadelphia (PA): SIAM; 2014.

31. , . Dolfin: automated finite element computing. ACM Trans Math Softw. 2010;37(2):417–44.

32. [[32]], , , . Unified framework for finite element assembly. Int J Comput Sci Eng. 2009;4(4):231–44.

33. Gaul A. Krypy – a python toolbox of iterative solvers for linear systems, commit: 23500bc9; 2014. Available from: https://github.com/andrenarchy/krypy

34. Heiland J. tdpbcvals – python module and test suite for convection-diffusion problems with time dependent Dirichlet boundary conditions, v2.0. 2015. Available from: https://gitlab.mpi-magdeburg.mpg.de/heiland/timedp-bcvals

35. , . Boundary concentrated finite element methods. SIAM J Numer Anal. 2003;41(1):1–36.

Competing Interests

The authors declare no competing interests.

Publishing Notes

© 2015 P. Benner and J. Heiland. This work has been published open access under Creative Commons Attribution License CC BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Conditions, terms of use and publishing policy can be found at www.scienceopen.com.