Time-dependent Dirichlet Conditions in Finite Element Discretizations

For the modelling and the 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 ﬁnite-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 conﬁrm that standard ﬁnite 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.


Introduction
In practical applications, see [20,21] for examples from 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 modelling of control setups, cf. [22,Ch. 1], timedependent 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., [6,Ch. 2] and [35,Ch. 2.3]. Another rather obvious obstacle is that a standard choice of trial and test functions formulations implies a certain smoothness of the boundary data which may be impractical [35,Ch. 2.3].
For a general overview of the functional analysis for parabolic systems with Dirichlet boundary control, we refer to [6,26]. 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 [7,18,29], time-dependent formulations are considered in [13,23]. The existence and the approximation of very weak solutions is well understood [7]. The key ingredient is the proper approximation of functions at the boundary via a projection [7,13,18].
An alternative approach of relaxing the boundary constraint via a penalization term in Robin boundary conditions has been investigated in [4,9].
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 [5] 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 forṁ v(t) = g(t, v, u) (1) or, in the linear case, in the forṁ 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 to 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 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 then 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 [16,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.
The paper is organized as follows. In Section 2 we state the type of problems that we will consider both in an abstract setting and after a spatial discretization. In Section 3, we consider approaches that reformulate the spatially discretized equations into the desired form. In Section 4, we discuss methods that reformulate the abstract equations such that a spatial discretization is a system of distributed type (1). In Section 5 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 modelling 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 Ω ∈ R d , d ∈ {2, 3}, be a bounded and regular domain such that the trace theorem as formulated in [8,Thm. 3.1] applies. Let Γ be its boundary. We define the Sobolev spaces V := W 1,2 (Ω) and H := L 2 (Ω) and the dual space V of V with respect to the continuous embedding of V in H to get We also introduce abbreviations for the trace spaces, cf.  Let be the trace operator as defined, e.g., in [12,Thm. 1.5].
We state the prototype of the continuous problem holds for almost all t ∈ (0, T ), and so that υ(0) = υ 0 in H.
The system of abstract equations (4) contains common weak formulations of PDEs that model physical phenomena, cf. [33]. We will not address time regularity here and, thus, leave the properties of the mappings t → A(t, υ(t)) and, e.g., t →υ(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 Ω ∈ R d , a diffusion parameter ν, a convection wind β, with β(x, t) ∈ R d for time t > 0 and x ∈ Ω, an initial value ρ 0 , and a function g, with g(t) : Γ → R prescribing the boundary conditions, find a function ρ of space and time that satisfiesρ and ρ(0) = ρ 0 .
In standard weak formulations, assuming υ ∈ V := W 1,2 (Ω), Problem 2.2 is of the type of Problem 2.1, with, e.g., A defined via for all φ ∈ V and with ∂ ∂n denoting the normal derivative.
Note that there are other possible choices for a weak formulation [23]. The boundary condition in Equation (4), viewed as a constraint, can also be incorporated using the dual operator of γ : V → Q and a so called Lagrange multiplier. Then, under certain smoothness and consistency conditions [14], Problem 2.1 is equivalent to Problem 2.3. Let T > 0 and consider A : (0, T ) × V → V . For F ∈ L 2 (0, T ; V ), υ 0 ∈ H, and U ∈ L 2 (0, T ; Q ), find υ with υ(t) ∈ V andυ(t) ∈ V and Λ with Λ(t) ∈ Q, a.e. on (0, T ), so thatυ holds for almost all t ∈ (0, T ), and so that υ(0) = υ 0 in H.
Note that, in general, the Lagrangian multiplier resides in the dual space of the constraint. In the considered case, where Q = [W 1 2 ,2 (Γ)] d is a Hilbert space, we can deliberately identify Q = Q.

Spatially Discretized Equations
We consider a generic spatial discretization of the introduced equations. Let V ⊂ V be a finite dimensional subspace spanned by the basis functions {ψ i } nv i=1 . 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 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 Γ ⊂ V is spanned by the basis functions {ψ i } nv i=n I +1 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 υ(t) ∈ V is to be approximated by a finite dimensional function v(t) ∈ V or the vector v(t) ∈ R nv containing the coefficients of the expansion in the considered basis. We will assume that v = (v I , v Γ ) is partitioned, with v I being associated with V I and v Γ 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 v ∈ V , we will identify v I and v Γ with their coefficient vectors of the expansion in (8) and simply write 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 := ψ i , ψ j H i,j=1··· ,nv with respect to the test space, and, once more, with respect to the trial space, where M II := ψ i , ψ j H i,j=1··· ,n I and M IΓ := ψ i , ψ j H i=1,··· ,n I ,j=n I +1,...nv 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 ) × R nv → R nv to A as and A I : (0, T ) × R nv → R n I as its restriction to the test functions of the inner nodes, where, again, we have associated a vector v ∈ R nv with a function in V via (8). If A is linear, then its approximation on A : (0, T ) × R nv → R nv can be assembled as a matrix-valued function via with the partitions A I , A II , and A IΓ as they were defined for M in (9).
The discrete approximation f : (0, T ) → R nv to the right-hand side F : (0, T ) → V is given as 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 and where u ∈ R nv−n I is a vector that will contain the current value of the boundary control at the given locations. As defined in (10), the operator G picks out the boundary dofs of a function v ∈ V and assigns the control values. Note, however, that there are other discrete approximations to the trace operator γ that consider, e.g. the inner product in Q or include suitable projections [7,18]. Thus, if we assume v(t) ∈ 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 v , n I ∈ N, and n d := n v − n I and consider A I : (0, T ) × R nv → R n I , G ∈ R n d ,nv , and M I ∈ R n I ,nv as defined in the beginning of Section 2.2. For f ∈ L 2 (0, T ; R n I ), α ∈ R nv , and u ∈ L 2 (0, T ; holds for almost all t ∈ (0, T ) and v(0) = α.
For the system of Problem 2.3 with the multiplier, a possible spatial discretization defines a differential equation considering also the boundary parts, cf. [2,14]. It generically takes the form Problem 2.5. Let T > 0, n v , n I ∈ N, and n d := n v − n I and consider A : (0, T ) × R nv → R nv , G ∈ R n d ,nv , and M ∈ R nv,nv as defined in the beginning of Section 2.2. For f ∈ L 2 (0, T ; R nv ), α ∈ R nv , and u ∈ L 2 (0, T ; R n d ) find v : (0, T ) → R nv and λ : (0, T ) → R n d , so that hold for almost all t ∈ (0, T ) and v(0) = α.
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 ∈ R n I ,nv and write (11) as 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 → H , with H = L 2 (Ω) where boundary conditions do not play a role.

Rewriting the Spatially Discretized Equations
In this section, we consider the spatially discretized equations introduced in Section 2.2. 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 [27], does not simply lead to a system of the form (1). Consider the linear formulation (13) of Problem the 2.4 with the assignment of the boundary conditions as in (13b): Then, with the partitioning of M I and A I as in (9), the state equation reads System (15) is not of the form (2) because of the appearance ofu.
Remark 3.1. One can define a new input asũ :=u and consider the system This approach uses a so called dynamical controller, that is defined via a differential relation. As pointed out in [3], 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 (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 τ , and subsequently discretizing in space, leads to systems of type where the superscript k relates to the values of the functions at the k-th time instance. If at the previous time instance v k Γ = u k , then the direct assignment at the current instance gives which can be seen as time-discrete approximation to (15).

Lifting of the Boundary Conditions
These approaches base on a liftingṽ that fulfills the boundary conditions for all time and the decoupling of the solution v = y +ṽ, see [32] for an example with linearized Navier-Stokes equations. We consider the linear time-invariant case (14) and assume that f = 0. At time t ∈ [0, T ], we define a lifting asṽ Then, considering (8) with v = y +ṽ and splitting A I and M I as in (9), we find that y Γ = 0 and we obtain the relation We use the abbreviationĀ II = M −1 II A II and, with the well-known solution representation, we obtain that After an integration by parts, we find that with Note, however, the dependency of the initial value on u(0) in (17). Remark 3.3. The dependency of the initial value on u(0) is due to the ansatz that assumes smoothness ofṽ, 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 which leads to M Iv (t) = 0 for all time t, the application for nonlinear systems is straight forward. Considering again, y = v −ṽ, and the nonlinear case of Problem 2.4, one arrives at the ODE Again, the actual solution is easily obtained by a backwards substitution v I = y I +ṽ I = y I − M −1 II M IΓ 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 ever 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: with the Lagrangian multiplier λ. 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 [15]. 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 = 0 I and M is symmetric positive definite. Thus, we can define With this, system (18) can be equivalently [16] andv Note that the differential equation (21) is of type (1).
With M P = P T M , in the linear case, we can write the differential equation for v i as In the nonlinear case, the input appears inside the nonlinearity.
Remark 3.7. Since n d n v , 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.  (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 αλ(t), 0 < α 1, to the left hand side of Equation (18b), one can solve for the multiplier and eliminate it from the differential part: This approach is known as penalty scheme and pressure penalization in the numerical integration of multibody and Navier-Stokes systems, respectively, cf., e.g., [10,31]. The method is straight forward to implement but comes with the need of a proper choice of the penalization parameter. The main difficulty is that a small parameter α 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 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 where λ : Q → V and α is a small penalization parameter. Then, for various choices of λ and Q, various weak incorporations of the Dirichlet conditions are realized. For example, defining λ through for a q ∈ Q and for any φ ∈ V, one obtains the penalized Robin approximation described in Section 4.3 below.

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 Φ = W 2,2 (Ω) ∩ W 1,2 0 (Ω) and consider the diffusion equation (5) with β = 0. We call v a very weak solution if for all φ ∈ Φ, cf. [18]. The abstract equations (23) indicate 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 which can be numerically approximated with standard finite element spaces of sufficiently high regularity.
With the nonconforming ansatz spaces V ⊂ W 1,2 0 (Ω), an approximation to (23) that is readily realizable, is given through v ∈ V that satisfies for all φ ∈ V , cf. [25,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 2 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) lead to decent approximations. The numerical approximation to ultra weak solutions of elliptic problems with boundary control as proposed in [7] uses a finite dimensional ansatz space V ⊂ H 1 (Ω) and as the test space W := V ∩ H 1 0 (Ω), see [18] and see [13] for the extension to parabolic problems. The elliptic case then reads, find v ∈ V , such that where Π V is the L 2 projection in L 2 (Γ) onto the grid induced by the triangulation that defines V . Using the spaces and formulations of (25) for a spatial discretization of a parabolic problem, one obtains a system that is the same as (14) apart from the appearance of the projector Π V in the boundary term (14b). Anyways, the elimination of the boundary nodes will lead to a system like (15) that includes the time derivativė u of the control. A discontinuous Galerkin ansatz for the time discretization as used in [13] includesu implicitly in the same way as Rothe's method, cf. Remark 3.2.
Remark 4.1. Since the known numerical approaches that base on (25) 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. Section 3.2, and the projection approach, cf. Section 3.3.1, readily apply to the formulation of the boundary term that includes the projector Π 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,13,18].

Nitsche Variational Formulation
A variant of the standard weak formulation of the pure diffusion, cf. (5) with β = 0, as proposed in [30] for the stationary Poisson equation reads for all φ ∈ Φ = W 1,2 (Ω). The formulation is derived by considering the cost functional with a parameter c γ and the first order optimality conditions for J (w − v) → min, where v is the solution of the stationary Poisson problem with nonhomogeneous Dirichlet boundary conditions. If for a given mesh c γ is chosen sufficiently large, namely c γ ≈ h −1 where h is a characteristic length of the triangulation, then the discretized optimization problem is convex [30, Eq. (12)]. For (26), a standard discrete formulation leads to an equation of type (2) with A and B explicitly given, see [34]. Cf. also [25,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 v ≈ α ∂v ∂n 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 α → 0 has been shown in several contexts, cf. [4] 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 modelling 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. Section 3.1 • lift -lifting of the boundary conditions via split mass matrix -cf. Section 3.2 • proj -incorporation of the constraint via Lagrange multiplier and projectionscf. Section 3.3.1 • pena -penalization of the constraint -cf. Section 3.3.2 • ncul -nonconforming approximation of ultra week solutions -cf. Section 4.1 • nits -approximation via the Nitsche variational-formulation -cf. Section 4.2 • pero -relaxation via Robin approximation -cf. Section 4.3 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. 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.

Test Setups
We consider several convection-diffusion setups on a two-dimensional square domain.
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
find approximations to the scalar function ρ satisfyinġ 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 2   The third test case is the same as the second but with an additional reaction source term r(ρ) = ρ(1 − ρ) in the dynamical equation. This source term r is positive for values of 0 ≤ ρ ≤ 1 and negative elsewhere. Thus, for values of ρ > 0 the reaction pushes ρ towards ρ = 1, cf. Figure 2 (b).
The considered system, for t ∈ (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 ρ satisfyinġ For all test cases, the spatial discretization is done on a uniform criss-cross triangulation described by the parameter N h = 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 τ ≈ 1 τ 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 P e := β(t) h/ν is smaller than 1. Thus, we can expect reliable approximations without additional, e.g. upwind, stabilization [24].
For the spatial discretization we use the Python interface dolfin [28] to the finite element software suite Fenics [1]. 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 implicit 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 errornorm 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 GMRES method using the implementation of the python module krypy [11]. The code used can be obtained from the author's public git repository [17].
By ρ pcg hN h ,τ Nτ , we denote the approximation to the solution of (28) with the discretization parameters N h , N τ , and cg. By e pcg hN h ,τ Nτ , we denote the approximation error e pcg hN h ,τ Nτ := ρ pcg hN h ,τ Nτ − ρ 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 τ = 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 τ . This justifies the choice of N τ := 240 as the reference discretization for further error comparisons.
The errors e pcg hN h ,τ 120 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 Section 5.5. For the linear elements, also ncul converges quadratically although with an error that is slightly larger than the one reported for proj and lift.
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 are the oscillation that occur 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 unevitable 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 γ = ν α 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 α.
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  Table 1: (Time space convergence of dias for linear elements, cf. Section 5.
2) The approximation error e pcg hN h ,τ Nτ with ρ ref =ρ p2 h96,τ 120 scaled by the inverse of e p1 h6,τ 30 = 1.119 · 10 −1 (left) and e p2 h6,τ 60 = 3.201 · 10 −2 (right) for varying space and time discretizations and for polynomial degree cg = 1 (left) and cg = 2 (right) for Test Case 1. Cf. also Figure 3(a,b) illustrating the convergence in space for the finest time discretization (the rightmost columns in the tables).  Table 2: (Time space convergence of dias for quadratic elements, cf. Section 5. 2) The approximation error e pcg hN h ,τ Nτ with ρ ref =ρ p2 h96,τ 120 scaled by the inverse of e p1 h6,τ 30 = 4.349·10 −4 (left) and e p2 h6,τ 60 = 8.551·10 −05 (right) for varying space and time discretizations and for polynomial degree cg = 1 (left) and cg = 2 (right) for Test Case 2. Cf. also Figure 3(c,d) illustrating the convergence in space for the finest time discretization (the rightmost columns in the tables). 3) The error e p1 hN h ,τ 120 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-h) to Test Case 3.
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 −4 .
The increase in the approximation error is mainly due to the increase of the magnitude of the right hand side that scales with 1 α . In fact, having solved the exemplary linear system Ax = f up to a relative residual of tol, one has that which means that for larger right hand sides f, the absolute residual Ax − f can be larger. A remedy is to control the absolute residual which can be done by correcting the provided relative residual by a factor tolcor = min{ 1 f , 1}, 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 1 α , one recovers the approximation properties of exact solves over the whole range of α, cf. Figure (6g-i) and (4d-f).

GMRES Performance
In this test setup, we investigated how the different but mainly equivalent formulations of the same problem affects the performance of an iterative solver. Therefore, we fixed the time and space discretization N h = 48 and N τ = 120 and, for polynomial degrees 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 Table 3 and Table 4.
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 base the two latest values. The parameter α was set to α = 1 for pena and α = 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 time-step 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 3. 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 Section 5, we have mentioned that the method of manufactured solutions is not suitable for boundary controlled processes. This is intuitively clear   since for ever 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 Case 1 and Test Case 2 (see Section 5.1) but with an additional volume force in (28a) corresponding to the constructed solution with u as in (27). The solution ρ ref is constructed such that at Γ 0 it coincides with the boundary control function defined in (27) and such that it is zero at the remaining boundaries. Also, it holds that ∂ρ ref ∂ν Γ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. Table 5, Table 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 3(b,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.

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 [19]. 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.