#### 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:

or, in the linear case, in the form: 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: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.

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 {\mathbb{R}}^{d},\text{\hspace{0.17em}}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:

Let

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* $\mathcal{F}\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*:

*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 {\mathbb{R}}^{d}$, *a diffusion parameter* $\nu $, *a convection wind β*, *with* $\beta \left(x,t\right)\in {\mathbb{R}}^{d}$ *for time t* > 0 *and* $x\in \Omega $, *an initial value ρ*_{0}, *and a function* $g$, *with* $g\left(t\right):\Gamma \to \mathbb{R}$ *prescribing the boundary conditions, find a function ρ* of *space and time that satisfies*:

*and ρ*(0) =

*ρ*

_{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:

*L*

^{2}(Γ), 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* $\mathcal{F}\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*:

*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\text{\hspace{0.17em}}=\text{\hspace{0.17em}}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:

where ${V}_{I}=\text{\hspace{0.17em}}\text{span}\text{\hspace{0.17em}}{\left\{{\psi}_{i}\right\}}_{i\text{\hspace{0.17em}}=\text{\hspace{0.17em}}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\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{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 {\mathbb{R}}^{{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:

*V*. If only Dirichlet conditions are posed, the standard test space is

*V*. Otherwise, all boundary dofs that are not fixed by a Dirichlet condition are included in the test space.

_{I}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:

Similarly, we define the discrete approximation $A:(0,T)\times {\mathbb{R}}^{{n}_{\upsilon}}\to {\mathbb{R}}^{{n}_{\upsilon}}$ to $\mathcal{A}$ as:

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

*A*,

_{I}*A*, and

_{II}*A*

_{I}_{Γ}as they were defined for

*in Equation (9).*

_{M}The discrete approximation $f:(0,T)\to {\mathbb{R}}^{{n}_{\upsilon}}$ to the right-hand side $\mathcal{F}:(0,T)\to \mathcal{V}\prime $ is given as:

*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:

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 {\mathbb{R}}^{{n}_{\upsilon}}\to {\mathbb{R}}^{{n}_{I}}$, $G\in {\mathbb{R}}^{{n}_{d},{n}_{\upsilon}}$, *and* ${M}_{I}\in {\mathbb{R}}^{{n}_{I},{n}_{v}}$ *as defined in the beginning of the section, Spatially discretized equations*. *For* $f\in {L}^{2}(0,T;{\mathbb{R}}^{{n}_{I}})$, $\alpha \in {\mathbb{R}}^{{n}_{\upsilon}}$, *and* $u\in {L}^{2}(0,T;{\mathbb{R}}^{{n}_{d}})$ *find* $\upsilon $ *with* $\upsilon (t):(0,T)\to {\mathbb{R}}^{{n}_{\upsilon}}$, *so that*:

*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 {\mathbb{R}}^{{n}_{\upsilon}}\to {\mathbb{R}}^{{n}_{\upsilon}}$, $G\in {\mathbb{R}}^{{n}_{d},{n}_{\upsilon}}$, *and* $M\in {\mathbb{R}}^{{n}_{\upsilon},{n}_{\upsilon}}$ *as defined in the beginning of the section Spatially discretized equations. For* $f\in {L}^{2}(0,T;{\mathbb{R}}^{{n}_{\upsilon}})$, $\alpha \in {\mathbb{R}}^{{n}_{\upsilon}}$, *and* $u\in {L}^{2}(0,T;{\mathbb{R}}^{{n}_{d}})$ *find* $\upsilon :(0,T)\to {\mathbb{R}}^{{n}_{\upsilon}}$ *and* $\lambda :(0,T)\to {\mathbb{R}}^{{n}_{d}}$, *so that*:

*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 {\mathbb{R}}^{{n}_{I},{n}_{\upsilon}}$ and write (4) as:

*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):

Then, with the partitioning of*M*and

_{I}*A*as in Equation (9), the state equation reads:

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

*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:

*k*th time instance. If at the previous time instance ${\upsilon}_{\Gamma}^{k}={u}^{k}$, then the direct assignment at the current instance gives: 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:

*A*and

_{I}*M*as in Equation (9), we find that ${y}_{\Gamma}=0$ and we obtain the relation:

_{I}After an integration by parts, we find that:

*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:

*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:

*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:

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:

With *MP* = *P ^{T}M*, in the linear case, we can write the differential equation for ${\upsilon}_{i}$ as:

*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:

*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):

*α*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:

*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:

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

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:

*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) 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:

*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:

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:

*α*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 {\mathbb{R}}^{2}$ be the computational domain with the spatial coordinates $x=({x}_{0},{x}_{1})$. Let Γ be the boundary with parts Γ_{0} to Γ_{3} 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.

The quantity *ρ* is seeded into the domain at Γ_{0}, where we enforce the time-dependent Dirichlet conditions:

_{1}, Γ

_{2}, and Γ

_{3}, 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:

*ρ*satisfying: 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.

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

*ρ*satisfying: 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:

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

*ρ*

_{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.

${\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 |

${\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.

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*.

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^{−4}.

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:

*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{\hspace{0.17em}}\text{=}\text{\hspace{0.17em}}\text{min}\left\{{\scriptscriptstyle \frac{\text{1}}{\Vert f\Vert}},\text{\hspace{0.17em}}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.

##### 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.

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.

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:

*ρ*

_{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.

${\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 |

${\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 |

#### 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.