A posteriori error estimation for parabolic problems with dynamic boundary conditions

This paper is concerned with adaptive mesh refinement strategies for the spatial discretization of parabolic problems with dynamic boundary conditions. This includes the characterization of inf-sup stable discretization schemes for a stationary model problem as a preliminary step. Based on an alternative formulation of the system as a partial differential-algebraic equation, we introduce a posteriori error estimators which allow local refinements as well as a special treatment of the boundary. We prove reliability and efficiency of the estimators and illustrate their performance in several numerical experiments.


Introduction
Within this paper, we consider a linear parabolic model problem with so-called dynamic boundary conditions in a bounded and polyhedral (Lipschitz) domain with initial condition u(0) = u 0 and right-hand sides f , ĝ. Therein, we write for the weighted Laplacians with uniformly positive diffusion coefficients α ∈ L ∞ (Ω) and κ ∈ L ∞ (Γ).
Condition (1b) itself is a differential equation and enables the reflection of effective properties on the boundary of the domain [20,26].To be more precise, this means that the momentum on the boundary is taken into account during the modeling process rather than being neglected as for standard boundary conditions.Because of this flexibility, this and similar models have attained increasing popularity nowadays, cf.[26,37,44].Although well-understood from a theoretical point (see, e.g., [22,42]), the numerical approximation is only tackled by a small number of papers; see [31,45] as well as the early work [21].
Within this paper, we benefit from a special formulation of (1) as partial differential-algebraic equation (PDAE); see [4,32] for an introduction.This means that we interpret the equations as a coupled system with an additional variable acting only on the boundary; see [34,Ch. 5.3] as well as [3].The resulting PDAE approach has already proven advantageous from a numerical point of view in the context of heterogeneous boundary conditions [6], phase separation models such as the Cahn-Hilliard equation [7], and the design of bulk-surface splitting schemes [5,8,16].The main instrument of the PDAE approach is the possibility to apply different discretization schemes in the bulk Ω and on the boundary Γ, e.g., based on different meshes T Ω and T Γ .This is also exploited in the present paper and is of special interest in applications where the solution oscillates rapidly on the boundary.Such oscillations also call for adaptive mesh refinements in order to improve the computational efficiency.Adaptive strategies, on the other hand, are based on a posteriori error estimates, which were first introduced for elliptic problems; see [2,36,43].
The paper is structured as follows.After the introduction of the PDAE model in Section 2.1, we shortly consider the temporal discretization and introduce a stationary model problem with a saddle point structure.Afterwards, we discuss inf-sup stable finite element schemes for this model problem in Section 2.3.
Based on classical finite element literature, we construct (local) a posteriori error estimators in Section 3.
Of special emphasis is the possibility to distinguish needed refinements of the bulk mesh on the boundary (denoted by T Ω | Γ ) and the boundary mesh T Γ .Moreover, we prove that these estimators are reliable as well as efficient, meaning that the estimators are indeed in the range of the actual error.We complete the paper with a number of numerical experiments in Section 4. This includes stationary problems on different polyhedral domains but also a dynamic problem such as (1).

Notation
Throughout the paper we write a ≲ b to indicate that there exists a generic constant C, independent of any spatial and temporal discretization parameters, such that a ≤ Cb.

Weak Formulation and Discretization
In this preliminary section, we introduce the weak operator formulation of the model problem (1) in PDAE form.Moreover, we discuss the discretization in time and space.

Weak formulation
With the introduction of an auxiliary variable p := tr u (with tr denoting the usual trace operator) as proposed in [3] and the spaces the weak form of (1) can be expressed as the following PDAE: Find u : This formulation includes the differential operators as well as the trace operator B u : V → M * and the canonical embedding We would like to emphasize that (2) depicts an extended formulation where the connection of u and its trace p is explicitly stated in the form of a constraint.In contrast to the classical formulation (1), the PDAE (2) includes a Lagrange multiplier as additional variable which satisfies λ = ∂ α,ν u if the solution is sufficiently smooth.Moreover, the system comes together with initial data u(0) = u 0 and p(0) = p 0 , which is called consistent if they coincide on the boundary, i.e., if B u u 0 = B p p 0 .The existence of a unique (distributional) solution of system (2) for consistent initial data follows from the inf-sup stability of the constraint operator together with the Gårding inequality of the differential operators; see [3] as well as [19].
At this point, we would like to emphasize that a direct spatial discretization of (2) leads to a differentialalgebraic equation (of differentiation index 2).In general, this implies numerical difficulties within the temporal discretization, cf.[27,28].Since the right-hand side of the constraint (2b) is homogeneous, however, these instabilities do not occur for consistent initial values and an index reduction is not needed, cf.[27, p. 33].
Within this paper, we follow the Rothe method for the discretization of system (2), i.e., we first discretize in time and then in space [30].The time discretization is shortly discussed in the following subsection, whereas we will focus on the spatial discretization using adaptive finite elements in each time step in Section 2.3.

Temporal discretization
We first consider the application of the implicit Euler scheme to (2) with constant step size τ.Given approximations u n ∈ V and p n ∈ Q of u and p at time point t n := nτ, respectively, we seek u n+1 ∈ V , p n+1 ∈ Q, and λ n+1 ∈ M as the solution of Therein, f (and similarly g) is some right-hand side which depends on the previous approximation u n (respectively p n ).
Remark 2.1.For semilinear applications, i.e., if we allow f and g to depend on u and p, respectively, we obtain the same linear structure as long as the nonlinear terms are treated explicitly.
In summary, an implicit Euler discretization leads to the following (stationary) problem, which needs to be solved in every time step: Given f ∈ V * and g ∈ Q * , find u ∈ V , p ∈ Q, and λ ∈ M such that This motivates the model problem, which we will consider in Section 3, namely Note that this model is written in its strong form because of which the Lagrange multiplier (being equal to the weighted normal trace) disappears in the first equation.Moreover, the newly introduced parameter σ > 0 corresponds to τ −1 in the case of an implicit Euler discretization.with an invertible coefficient matrix A ∈ R s,s , which defines an algebraically stable Runge-Kutta scheme.
According to [46,Sect. 8.3 f.], this leads to a stable approximation.Furthermore, we assume that −A is globally stable, which is equivalent to the existence of a symmetric matrix D ∈ R s,s such that D and DA −1 + A −T D are positive definite [33,Th. 13.1.1].Then, the calculation of the internal stages u n+1 ∈ V s , p n+1 ∈ Q s , λ λ λ n+1 ∈ M s can be treated similarly to the implicit Euler method considered above.In comparison to (3), u n+1 (and analogously p n+1 and λ n+1 ) needs to be replaced by u n+1 .The operators are replaced by their Kronecker product with D, e.g., K turns into D ⊗ K .Moreover, the operator DA −1 ⊗ id pops up in front of u n+1 and p n+1 .Finally, f (and analogously g) contains the approximation of the previous time point u For the translation of the upcoming analysis given in Section 3 to the Runge-Kutta case, one uses that the assumptions on the matrix D imply Examples of diagonal matrices D for the Gauss-Legendre, Radau IA, and Radau IIA methods can be found at [28, p. 220 f.].

Stable spatial discretization
We first collect a number of standard finite element spaces, which will be used in the following.Given a regular triangulation T Ω of the computational domain Ω, which is assumed to be polyhedral, the space of (discontinuous) piecewise polynomials of degree k ≥ 0 is denoted by The subspace of functions which are additionally elements of H 1 (Ω), is denoted by P k (T Ω ).It is wellknown that this space can be characterized by see, e.g., [10,Th. II.5.2].Given a regular triangulation T Γ of the boundary Γ we denote analogously the space of (possibly discontinuous) piecewise polynomial functions and its subspace of continuous functions by P d k (T Γ ) and P k (T Γ ), respectively.For the construction of stable finite element schemes, we further introduce the space of bubble functions.In two space dimensions, the edge-bubble function ψ E equals the scaled product of the two nodal basis functions corresponding to the nodes of the edge E. For d = 3, ψ E denotes the face-bubble (E being a face in T Ω ) and equals the scaled product of the three corresponding nodal basis functions, cf.[43].This leads to the space ψ E is an edge/face-bubble for E ⊆ Γ ⊆ P ℓ+d (T Ω ).Now let V h , Q h , and M h denote finite-dimensional subspaces of V , Q, and M .Due to the saddle point structure of the model problem, these spaces need to satisfy a discrete inf-sup condition, which reads inf for some positive constant β being independent of the mesh sizes.Here, ⟨ • , • ⟩ Γ denotes the dual product in M , which only acts on Γ.In a first result, we show that this stability condition is independent of the choice of the space Q h .
Theorem 2.3 (Equivalence of inf-sup stable discretizations).Consider a conforming Galerkin scheme , which includes the approximation property of the family V h , i.e., inf v h ∈V h ∥u − v h ∥ V → 0 as h → 0 + for all u ∈ V .Then, the discrete inf-sup condition (5) is satisfied if and only if there exists a positive constant β (independent of the mesh size) such that inf Proof.(⇒) We prove the if-part by contrapositive and assume that ( 6) is not satisfied.Then, there exists Since the sequence is uniformly bounded, there exists a subsequence implies the existence of a u ∈ V with B u u = g.Due to the assumed approximation property of V h , there exists a sequence {u h } h ⊆ V with u h ∈ V h , u h ̸ = 0, and lim h→0 + u h = u in V .By the strong convergence of {u h } h and the weak convergence of {λ h ′ } h ′ , we have Hence, λ * = 0 holds and the entire sequence {λ * h } vanishes weakly in M = H − 1 /2 (Γ), cf.[24, Ch.I, Lem.5.4].Moreover, the sequence vanishes strongly in Q * = H −1 (Γ), since the weak limit is unique and the embedding where we have used that B p is the embedding operator of . Thus, the inf-sup expression cannot be bounded uniformly from below by a β > 0.
(⇐) The only-if-part follows immediately by sup Remark 2.4.In the limit case κ ≡ 0, there is no differential operator acting on the boundary and the boundary conditions are called locally reacting.The corresponding solution spaces then read In this setting, one can show analogously to Theorem 2.3 that conforming finite element spaces are inf-sup stable if and only if the subproblem with A direct consequence of Theorem 2.3 is that it is sufficient to consider the inf-sup condition with p h = 0. Hence, one example class of inf-sup stable discretizations reads We would like to emphasize that the mesh used for the definition of M h may be different to T Ω | Γ as long as it is synchronized with the mesh used for the bubble functions.Moreover, the scheme is inf-sup stable; see [6,Prop. 3.4].The corresponding generalization to higher polynomial degrees reads as follows.
Lemma 2.5.The conforming finite element spaces satisfy the discrete inf-sup condition (5) for arbitrary k ≥ 1 and Q h .
Proof.The case k = 1 is shown in [6,Prop. 3.4].The statement for general k can be proven equivalently.
Here, we use that there exists an extension operator from Recall that all presented stable schemes still have the freedom to set the discrete space Q h .In the following, Q h will be a standard finite element space but defined on another triangulation T Γ .Having in mind applications with strong fluctuations or heterogeneities on the boundary, T Γ usually equals a refinement of T Ω | Γ .
Remark 2.6.The choice V h = P 2 (T Ω ), Q h = P 1 (T Γ ), and M h = P 1 (T Ω | Γ ) is stable and yields a priori rates of order h 2 for u and h for p and λ (in the respective H 1 -norms).Due to the different spatial dimensions, this then implies a balanced order of one in terms of the degrees of freedom.

A Posteriori Error Estimation
After the discussion on the temporal discretization and inf-sup stable finite element schemes, we now turn to the construction of efficient and reliable a posteriori error estimators.Following the Rothe method, the aim is an adaptive spatial discretization in each time step.Note that the formulation as coupled system in (2) amplifies the use of adaptive schemes as we can optimize T Ω (for the spaces V h , M h ) as well as T Γ (for Q h ).Throughout this paper, we assume T Ω and T Γ to be shape regular; see [10, Ch.II.5].Note that this automatically implies the shape regularity of T Ω | Γ .Within this section, we consider the two-dimensional setting but comment on necessary adjustments in the three-dimensional case.
By E Ω we denote the set of edges corresponding to the triangulation T Ω , which can be decomposed into E in Ω (interior edges) and E bd Ω (boundary edges).Moreover, h T and h E denote the mesh sizes for elements and edges of a triangulation T Ω , respectively.For a partition T Γ of the boundary, we write h I for the mesh size of these elements.For the upcoming analysis, we consider the following assumption.(i) The discretization scheme is inf-sup stable and it holds that P 1 (T Ω ) ⊆ V h , P 1 (T Γ ) ⊆ Q h , and Moreover, there exists a positive constant ρ > 0 such that h E ≤ ρh I holds uniformly for the mesh sizes h I and h E of every As explained in Section 2.2, we deal with the model problem (4).Working again with the weak form, we consider the system for test functions v ∈ V , q ∈ Q, µ ∈ M and the parameter σ which is related to the time step size.
Remark 3.2.The upcoming analysis can be performed similarly for semilinear applications; see [43,Sect. 5.2.3].We, however, focus on the linear case, since our main interest is the treatment of the two different meshes on the boundary.
To shorten notation, we define the product space Moreover, we introduce the bilinear form B : X × X → R, which is directly related to the model problem, by Note that we use the symbolic integral notation for the boundary terms as they equal the L 2 -inner product if the arguments are sufficiently smooth.Obviously, the bilinear form B is continuous.Furthermore, following the lines of [43,Prop. 4.68], one can show that B is inf-sup stable, i.e., there exists a positive constant With B at hand, the introduced weak form of the model problem can be written as where the At this point, we would like to recall that we are especially interested in cases where strong fluctuations occur on the boundary.This means, in particular, that we expect Q h to be a refinement of V h restricted to the boundary.Moreover, this motivates the use of adaptive methods, which are based on a posteriori error estimates.These are topic of the upcoming subsection.

Construction of local estimators
Following the methodology of [36,43], we aim to construct a posteriori error estimators which can distinguish necessary refinements for the approximations of u| Γ and p, respectively.Within this section, the triple [u, p, λ ] ∈ X denotes the weak solution to the model problem ( 4) and [u h , p h , λ h ] ∈ X h its discrete counterpart.As a consequence of (7), we have the error bound In order to bound the right-hand side, we integrate by parts and obtain Recall that E in Ω and E bd Ω denote the set of interior and boundary edges of the triangulation T Ω , respectively.Moreover, N Γ equals the set of nodes of T Γ and [ • ] E the jump along an edge E.

Making use of the mentioned Galerkin orthogonality, namely
we can add arbitrary discrete test functions.Hence, including v h ∈ V h as any quasi-interpolant of v (such as the Clément interpolant [15]) and q h ∈ Q h as the (pointwise) interpolant of q, we can estimate Here, we have used standard (quasi) interpolation results as presented, e.g., in [13].This includes the well-known element and edge patches, denoted by ω T and ω E , respectively.Further note that the difference tr u h − p h occurs globally, since the H 1 /2 (Γ)-norm is not additive.For a local version, we use that T Γ is a refinement of T Ω | Γ and restrict ourselves to a bisection strategy for the refinement of T Γ .This leads to a local estimate which only depends on the shape regularity of the bulk mesh T Ω .Before we show the local estimate, we fix the refinement strategy in the following assumption.Assumption 3.3.Given the restriction of the bulk mesh to the boundary T Ω | Γ , there exists a sequence of refinements using bisection to obtain the boundary mesh T Γ .
We would like to emphasize that Assumption 3.3 is satisfied, if the initial meshes of Γ and Ω are refined by bisection and newest vertex bisection (NVB) [38], respectively.Lemma 3.4.Given Assumptions 3.1 and 3.3, it holds that with a hidden constant only depending on the polygonal domain Ω with its boundary Γ and the shape regularity of the triangulation T Ω .
Proof.By the assumed bisection strategy, there exists a refinement T Ω of T Ω such that T Γ = T Ω | Γ with the shape regularity of T Ω only depending on the shape regularity of T Ω ; see [ Here, h Γ,max (respectively h Γ,min ) denotes the maximal (minimal) mesh width of the triangulation T Γ .Lemma 3.4 states that a refinement based on Assumption 3.3 yields a quasi-uniform triangulation for which this ratio is bounded.
Finally, we define the following local quantities, which will serve as local error estimators, namely The remainder of this section is devoted to the verification of the reliability and efficiency of these estimators.
Remark 3.7.In the above derivation of the error estimators, we have used that functions in H 1 (Γ) are continuous for a one-dimensional boundary Γ of a two-dimensional domain Ω.This allows us to neglect the jumps of p over the nodes of T Γ .If Ω is a subset of R 3 , however, we have to consider these jumps, which are now jumps over edges.As a consequence, the estimators in (9) need to be extended by for D ∈ E Γ (10) with E Γ being the set of edges of T Γ and h D the length of an edge D ∈ E Γ .

Reliability and efficiency
The reliability of the introduced a posteriori estimators follows directly by construction and is summarized in the following theorem.
Theorem 3.8 (Reliability of the error estimator).Let η T , η E , and η I be defined as in (9) and let Assumptions 3.1 and 3.3 be satisfied.Then the discretization error can be bounded by where the hidden constant only depends on the polynomial degrees used for the finite element spaces, the polygonal domain Ω with its boundary Γ, and the shape regularity of the triangulations T Ω and T Γ , and the coefficients α, κ, and σ .
Remark 3.9.The reliability in the three-dimensional case follows similarly to Theorem 3.8, where we have the additional term ∑ D∈E Γ η 2 D on the right-hand side, cf.Remark 3.7.
Theorem 3.8 guarantees that the error estimators provide an upper bound of the error.It remains to prove the efficiency of the estimators, i.e., the guarantee that they are (up to a constant) also a lower bound.
In the following lemmas, we always assume that Assumption 3.1 is satisfied.Similar to Theorem 3.8, all hidden constants will only depend on the dimension, the polynomial degrees of the finite element spaces, the polygonal domain with its boundary, the shape regularity of the meshes, the diffusion coefficients, and the constant ρ.We first consider the estimators η T and η I .
Lemma 3.10.The estimators η T and η I satisfy as well as These estimates further imply Proof.Estimate (11a) follows by the standard bubble function technique, see, e.g., [36,Lem. 2.6.10(1)].For (11b), we note that every segment . With this, one proves (11b) analogously to (11a), where one uses that for every ν ∈ H 1 0 (I) -extended by zero outside of I -it holds the interpolation inequality To obtain an upper bound on η E , we distinguish interior and boundary edges.
Lemma 3.11.For an interior edge E ∈ E in Ω it holds that Proof.The proof follows the lines of [2,Lem. 3.6].
Lemma 3.12.For a boundary edge E ∈ E bd Ω we define the boundary edge patch ω E,Γ by and T E as the (unique) element in T Ω with edge E.Then, it holds that and, hence, Proof.Combining the arguments of [36, Lem.2.6.12],[2,Lem. 3.6], and the one used for the derivation of (11b), we have and, hence, for the entire boundary For an estimate of the second summand of η E , we set µ with the characteristic function χ I .We now distinguish the two cases Then, the continuous and discrete constraint imply that By [1, Prop.2.2 & 4.8], there exists a constant only depending on the shape regularity of Combining the latter two estimates, we conclude that With Lemma 3.10 and appropriate Sobolev embeddings, we get the stated estimate for η E .Finally, the estimate of the sum follows by where the integrals over Γ are defined via local charts.
) be the (weighted) Clément operator Since µ vanishes outside of E, its quasi-interpolation P E µ is zero on the relative boundary ∂ ω E,Γ .Hence, P E µ can be extended by zero to a function in M h .Similar to the previous case, we obtain the estimate By [11,Lem. 3.1], we have ∥µ − P E µ∥ L 2 (ω E,Γ ) ≲ ∥µ∥ L 2 (E) as well as with constants only depending on the regularity of T Ω | Γ .Here, we have used that P E is symmetric with respect to the L 2 (ω The remainder of the proof follows the lines of the first case.
It remains to summarize the previous estimates.The combination of Lemmas 3.10, 3.11, and 3.12 yield the following efficiency result, showing that the error estimators are bounded from above by the actual error plus oscillation terms of the right-hand sides.
Theorem 3.13 (Efficiency of the error estimator).Given Assumption 3.1, the error estimators defined in (9) satisfy The hidden constant only depends on the polynomial degrees of the finite element spaces, the polygonal domain Ω with its boundary Γ, the shape regularity of the triangulations T Ω and T Γ , the coefficients α and κ, as well as the constants σ and ρ.
Remark 3.14.Theorem 3.13 shows that the estimated errors are bounded by the actual error plus data oscillation terms.Given a discretization satisfying Assumption 3.1 and right-hand sides in the respective H 1 -spaces, theses terms can be estimated by as shown in [10, Ch.II, Cor.7.7].Thus, in comparison with the error ∥[u − u h , p − p h , λ − λ h ]∥, the data oscillation terms are of higher order and may be neglected.
Remark 3.15.Analogously to η E in Lemma 3.11, the additional error estimator η D defined in (10) is bounded by for a three-dimensional polyhedral domain Ω.Here, ω D is defined with respect to T Γ as ω E to T Ω .Furthermore, Lemmas 3.10, 3.11, and 3.12 are also valid for Ω ⊆ R 3 , implicating again the efficiency of the error estimators.

Refinement and coarsening
With the shown reliability and efficiency of the error estimators (9), we can now discuss sound strategies for the generation of spatial meshes for the semi-discrete system (3).
For elliptic problems, it is well-known that a refinement using NVB based on a Dörfler marking [18] yields an optimal strategy [14].Hence, we apply this refinement strategy in every time step.However, since NVB is based on element marking rather than marking elements and edges, we attribute the estimator η E of an edge E among all η T and η I with T ⊆ ω E and I ⊆ E, respectively.The resulting error estimators are denoted by η T and η I , respectively.Then we use the Dörfler marking with a parameter θ ∈ (0, 1), i.e., we mark all elements in (minimal) sets M Ω ⊆ T Ω and M Γ ⊆ T Γ such that and refine every marked element using bisection or NVB, respectively.
For time-dependent problems, the region of the bulk or boundary with the largest estimated errors can change in every time step.Hence, taking the final meshes of the previous time step as the initial meshes for the current time step may lead to a fast growing number of degrees of freedom.In particular, a strongly refined discretization in a certain region may become obsolete after some time.Therefore, it makes sense not only to refine the meshes over time but to allow coarsening as well if needed; see [9,Sec. 4.4.4].Here, we follow the heuristically motivated coarsening strategy of [23,Sec. 7.2].For this, we mark elements T ∈ T Ω and I ∈ T Γ for a coarsening step if they satisfy for a given tolerance tol and a parameter γ ∈ (0, 1).Note that in the case of a semi-discretized system, the right-hand sides f n and g n include the approximations of the previous time step u n−1 h and p n−1 h , respectively.After a possible coarsening, the meshes at the current time point t n may not be refinements of the meshes from time t n−1 anymore.In order to represent u n−1 h and p n−1 h in these meshes, we consider their L 2 -best approximations.This produces additional errors in the right-hand sides, which can be estimated similarly to the oscillation data discussed in Remark 3.14.Hence, we obtain additional error terms of the form . To summarize, in every time step one first refines the bulk and the boundary mesh using NVB and bisection, respectively, in combination with the Dörfler marking strategy until the estimated error is below a given tolerance.This yields approximations of u, p, and λ for the current time point.Afterwards, we coarsen the meshes by repetitively marking elements T and I satisfying (13) until either no marked elements are left or only elements are marked where a coarsening would contradict Assumption 3.1.(ii).The resulting meshes are then used as starting point for the upcoming time step.

Numerical Experiments
In this final section, we demonstrate the performance of the newly introduced a posteriori error estimators applied to the stationary model problem (4) as well as a parabolic problem with dynamic boundary conditions.
In all examples, we consider a two-dimensional domain and diffusion coefficients α ≡ 1, κ ≡ 1.As refinement strategy, we use bisection as refinement for T Γ and NVB for T Ω as discussed in Section 3.3 with the implementation taken from [23].For the computation of the errors, we use a reference solution, which is computed on a refined mesh.More precisely, we consider in each step the current mesh obtained by the adaptive process with two additional uniform refinements.

Stationary problem on the unit square
In this first example, we consider the stationary model problem (4) on the unit square, i.e., on Ω = (0, 1) 2 .The right-hand sides are given by f ≡ 0.04, g(x, y) = xy cos(10πx) cos(10πy) and we set σ = 1 as well as θ = 0.75 for the marking process.The resulting convergence history for the scheme is illustrated in Figure 1.Note that the x-axis includes the sum of the degrees of freedom for u, p, and λ .The plot clearly shows the improvement caused by the adaptive process as we reach a convergence rate of 0.65 rather than 0.5 for a uniform refinement.Note that a rate larger than 0.5 is only possible due to the mixture of different schemes (and dimensions).Considering the single components, one observes for u, p, and λ convergence rates of 0.5, 1.0, and 1.5, respectively.Moreover, one can observe that error and estimator, which are defined as error are indeed of the same order as expected due to Theorems 3.8 and 3.13.
The resulting mesh of the adaptive process is shown in Figure 2. It clearly shows a refinement of the upper right corner for T Ω as well as for T Γ .This is due to the fact that the right-hand side g as well as the solution oscillate strongly in this region.

Stationary problem on the L-shape
The second example works on the L-shape, i.e., on a non-convex computational domain.Here, we consider the right-hand sides and again σ = 1, θ = 0.75.We compare two different stable finite element schemes.
First, we use piecewise linear elements as in (14).As shown in Figure 3, this yields the rate 0.375 for a uniform refinement and the optimal rate 0.5 in the adaptive case.Second, we consider a piecewise quadratic approximation for u and p in combination with a piecewise constant approximation of the Lagrange multiplier, i.e., Note that this is indeed a stable scheme according to the derivations of Section 2.3.The convergence results are shown in Figure 4 and demonstrate an even higher numerical gain of the adaptive procedure.
If the discrete multiplier space M h is set to P 1 (T Ω | Γ ) or P 2 (T Ω | Γ ), we observe the same rates.Hence, we omit the corresponding error plots here.
An illustration of the adaptive mesh refinement is given in Figure 5.Note that both meshes show refinements in the same regions.Nevertheless, u restricted to the boundary has 190 degrees of freedom, whereas p is defined on a mesh with 1030 degrees of freedom.This additional refinement of the boundary is only possible because of the special formulation of the system equations as coupled system.Without the separation of u| Γ and p, we would need a refinement of T Ω | Γ in order to get the same accuracy.This, however, would call for an significant increase of the degrees of freedom due to the higher topological dimension of the bulk.with homogeneous initial data.For the computation, we set [1,10] as time horizon and apply the implicit Euler scheme with constant step size τ = 1.5•10 −2 for the time stepping.For the spatial discretization, we consider piecewise quadratic elements as in (15).At every time point, we check if the estimated spatial error is smaller than 10 −6 .If not, we use the mentioned Dörfler marking strategy with parameter θ = 0.75 and refine the meshes until this condition is satisfied.Hence, no coarsening is used in this experiment.

Parabolic problem with dynamic boundary conditions
In the beginning, we start with 41, 16, and 8 degrees of freedom for u, p, and λ , respectively.At the end of the simulation, the numbers are 4912, 20674, and 301.The detailed development of the degrees of    As a proof of concept, we repeat this experiment with a coarsening strategy for the boundary mesh.As marking condition, we use η I ≤ 1/(4 • 10 6 |T Γ |); cf.Section 3.3.The associated degrees of freedom over time are shown in Figure 6 (right).Obviously, the numbers are always smaller compared to the simulation without coarsening and are minimal for t ≈ 1.5, 2.5, . .., i.e., whenever g(t) = 0.This coincides with the times where the degrees of freedom of p stay almost constant in the simulation without coarsening.At these points, a sound representation of p does not need a fine mesh and the coarsening works as intended.
We would like to mention that the simulation with coarsening takes much longer than the one without.The degrees of freedom in the present example, however, are still quite low and we expect a dominant performance if regions needing a very fine discretization change more frequently.
Finally, we present the approximated solution u at the final time point t = 10 in Figure 7.One can see that the solution is oscillating strongly at the boundary parts (0, 1) × {1} and {1} × (0, 1).To capture this behavior numerically, we need more degrees of freedom for p.On the other hand, the approximation of u gets along on a coarser mesh, since the solution is not so oscillatory in the bulk.

Conclusion
In this paper, we have introduced an adaptive procedure for the simulation of parabolic problems with dynamic boundary conditions.Based on a PDAE formulation, we are able to construct local error estimators, which distinguish necessary refinements of the bulk mesh along the boundary, i.e., of T Ω | Γ , and the boundary mesh T Γ .Moreover, we have characterized inf-sup stable finite element schemes for the spatial discretization.Within a number of numerical experiments, we have shown that the proposed algorithm reaches optimal convergence rates.In addition, one can observe that the boundary mesh gets more refined than T Ω | Γ , which yields notable computational savings.

Figure 1 .
Figure 1.Convergence history for the stationary model problem of Section 4.1.The dashed line in gray indicates order 0.5, whereas the solid gray line indicates order 0.65.

Figure 2 .Figure 3 .
Figure 2. Illustration of the bulk mesh T Ω (left) and the boundary mesh T Γ (right) after 40 adaptive refinement steps.

Figure 4 .
Figure 4. Convergence history for the stationary model problem of Section 4.2 using piecewise quadratic elements for u, p and piecewise constants for λ .The gray lines indicate orders 0.375 (dashed) and 1.0 (solid).

Figure 5 .
Figure 5. Illustration of the bulk mesh T Ω (left) and the boundary mesh T Γ (right) after 40 adaptive refinement steps.

Figure 6 .
Figure 6.Evolution of the degrees of freedom for u, p, and λ over time without (left) and with (right) coarsening of Q h .

Figure 7 .
Figure 7. Numerical solution u at the final time point T = 10 for the parabolic problem of Section 4.3.
Standard arguments such as in [12, Sect.II.2.2] imply the a priori estimate