{'fraction_non_alphanumeric': 0.03605121293800539, 'fraction_numerical': 0.012702156334231806, 'mean_word_length': 4.188079007166579, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Existence, uniqueness, and approximation of a fictitious domain formulation for fluid-structure interactions

In this paper we describe a computational model for the simulation of fluid-structure interaction problems based on a fictitious domain approach. We summarize the results presented over the last years when our research evolved from the Finite Element Immersed Boundary Method (FE-IBM) to the actual Finite Element Distributed Lagrange Multiplier method (FE-DLM). We recall the well-posedness of our formulation at the continuous level in a simplified setting. We describe various time semi-discretizations that provide unconditionally stable schemes. Finally we report the stability analysis for the finite element space discretization where some improvements and generalizations of the previous results are obtained.


Introduction
In this paper we summarize in a unified setting some results of our research on the modeling and the approximation of fluid-structure interaction problems. Our aim is to describe the dynamics of a solid elastic body immersed in a Newtonian incompressible fluid. Here, we consider the so called zero-codimension case, that is the solid and the fluid are both two-or three-dimensional. From the mathematical point of view, the interaction is described by different partial differential equations in the regions occupied by the fluid and the solid, coupled with suitable transmission conditions along the interface between the two. It is well known that the numerical approximation of fluid-structure interaction problems is challenging for several reasons: first of all the numerical method must track the movement of the structure and the corresponding computational grids should allow the evaluation of quantities defined on moving domains. In this context, the use of a Lagrangian framework is more suited for the simulation of the structure deformation, while the approximation of the fluid velocity and pressure is better performed by an Eulerian approach.
Another crucial issue related to the approximation of fluid-structure interactions is how to deal with the coupling of the two underlying models: monolithic approaches perform the simultaneous computation of the fluid and structure unknowns, while partitioned schemes combine different solvers in the two subregions with an iterative procedure. In general monolithic schemes require implicit nonlinear solvers and a careful trade off between superior stability properties and more demanding computational load.
The research in this framework is very active and is based on a wide literature, ranging from boundary fitted approaches which, typically, use the so called Arbitrary Lagrangian Eulerian method [38,28,40,29] to non fitted approaches which include, for instance, level set methods [21] Nitsche and XFEM methods [18,1].
Our model belongs to the latter family originating from the Immersed Boundary Method (IBM) [45,7] and evolved towards a fictitious domain approach in the spirit of [35,34,31,33,32,52]. Obviously, no method is the optimal choice for all cases and, depending on the particular situation, it could be preferable to make use of different approaches; our formulation has the advantage to be unconditionally stable in time [6,15] without the need of using fully implicit time schemes and, being based on non fitted meshes, can accommodate larger displacements. On the other hand, the coupling between fluid and structure models requires the evaluation of integrals that combine basis functions defined on different meshes. A solid mathematical analysis has been performed; we shall review some of the results in the following sections giving reference to the original papers when appropriate. Moreover, we extend the discretization of our model, allowing for more general choices of finite element spaces. We describe an incompressible solid immersed in an incompressible fluid; more general situation could be considered, involving compressible solids [13].
In Section 2 we recall the problem we are interested in, and introduce our fictitious domain formulation. Next, we analyze the continuous problem in Section 3 in a linearized setting, assuming that the motion of the solid is prescribed. Section 4 deals with the time discretization; the main result of this section is the unconditional stability of the evolution scheme. The space discretization is considered in Section 5 where a stability analysis is presented which leads to optimal convergence estimates for the steady state solution. Finally, Section 6 reports on several numerical tests that confirm the good behavior of our approach.

Model problem and fictitious domain formulation
The problem we want to address is easily explained in the following simplified setting. We consider a solid immersed in a fluid in two or three dimensions. At time t the solid is located in the region Ω s t ⊂ R d (d = 2, 3) which is the image of a reference configuration B through a mapping X : B → R d . The fluid occupies the region Ω f t ⊂ R d so that we are interested in a dynamic occurring in the union of Ω s t and Ω f t . A typical assumption is that, denoting by Ω the interior of the union of the closures of Ω s t and Ω f t , then Ω does not depend on t. This assumption is reasonable for several applications; in general Ω can be thought as a container where the dynamics takes place: for instance, the solid can be inside the fluid and far away from the exterior boundary of it, or the solid can touch one fixed part of the container. In this paper we deal with the first situation. We denote by Γ t the interface between fluid and solid, which can be defined as the interior of the intersection of Ω s t and Ω f t . The system is described by the fluid velocity u f and pressure p f , and by the solid position X. The velocity and the pressure depend on time and on the space Eulerian variable x ∈ Ω f t , while the position X depends on time and on the Lagrangian variable s ∈ B. In the fixed domain Ω we are using the Eulerian framework and the corresponding variable x. A point x of the domain Ω s t can be expressed at time t in the Lagrangian setting as The kinematic condition is expressed by the following relationship between the material velocity u s and X: where x = X(s, t). The deformation gradient is given by We denote by |F| its determinant. We consider an incompressible solid, so that |F| is constant in time; in particular, in the case when B is the initial configuration Ω s 0 of Ω s t , we have |F| = 1. In the incompressible fluid the Navier-Stokes equations describe the dynamics as follows (1) ρ where ρ f is the fluid density and σ f is the Cauchy stress tensor that reads ν f > 0 being the viscosity of the fluid and ε the symmetric gradient. We assume an incompressible viscoelastic material that can be described by a Cauchy stress tensor composed of two parts σ s = σ s f +σ s s : the first one is analogous to the fluid stress with the introduction of an artificial pressure p s , which is the Lagrange multiplier associated with the incompressibility, σ s f = −p s I + ν s ε(u s ), ν s > 0 being the body viscosity; the second term is related to the Piola-Kirchhoff elasticity stress tensor P via the Piola transformation σ s s = |F| −1 PF . The elastic part of the stress can be modeled using a potential energy density W (F, s, t) so that Taking all this into account, the equations describing the solid are in Ω s t , where ρ s is the solid density. The description of the model requires suitable transmission conditions enforcing the appropriate continuities of the velocity and of the Cauchy stress across the interface Γ t which can be stated as follows where n f and n s stand for the outward unit normal vectors to Ω f t and Ω s t , respectively. In conclusion, the system is described by (1), (2), (3), and the following initial and boundary conditions Before describing our variational formulation we recall some standard notation that we are going to adopt [41]. Given a domain D, the space D(D) is the space of infinitely differentiable functions with compact support in D, L 2 (D) is the space of square integrable functions on D, the standard Sobolev spaces are denoted by W s,p (D), where s ∈ R refers to the differentiability and p ∈ [1, +∞] to the integrability exponent. As usual, when p = 2 we use the notation H s (D). The corresponding norm is indicated by · s,D and the scalar product in L 2 (D) by (·, ·) D ; when no confusion arises we omit the indication of the domain D. In particular we will usually omit Ω, while we will indicate explicitly when quantities are defined on the domain B. L 2 0 (D) stands for the subspace of zero mean valued functions and H 1 0 (D) is the subset of functions in H 1 (D) with zero trace on ∂D. Given Banach spaces X and Y , the notation Y (0, T ; X) contains space-time functions that for almost all t ∈ 0, T are in X and that are in Y as functions from (0, T ) to X. Functional spaces of vector valued functions are indicated with boldface letters.
The main idea behind the fictitious domain approach that we are going to adopt, consists in extending the fluid variables inside the solid domain so that all involved quantities are defined in Ω (Eulerian variables) or B (Lagrangian variables). We started considering a fictitious domain model for a simplified interface problem [2,14] which has been extended to fluid-structure interactions in [6].
We denote by u and p the velocity and pressure in Ω, with the understanding that their restrictions to the two subdomains Ω f t and Ω s t coincide with u f , p f and u s , p s , respectively. With the aim of presenting a variational formulation of our problem, the condition u| Ω s t = u s will be enforced with the help of a bilinear form. Let Λ be a Hilbert space and c : Λ × H 1 (B) → R a continuous bilinear form with the property The variational formulation is described by making use of the following notation.
In our setting two possible definitions of c have been discussed in [6,9]: a natural choice is to consider as Λ the dual space of H 1 (B) so that c can be taken as the duality pairing that certainly satisfies the required properties; a second equivalent choice stems from interpreting the duality pairing as the scalar product in H 1 (B) by the Riesz representation theorem so that Λ = H 1 (B). More in detail, we have the following definitions 1. Λ 1 = H 1 (B) and c 1 : While the two definitions are equivalent for the continuous problem, they give rise to different discretizations. In the sequel we are going to use the generic notation Λ and c, while indicating explicitly one of the two cases when needed.
An analogous formulation, which is outside the topics of the present work, can also be used in the case of codimension one structures. We refer the interested reader to [6,9].
We end this section by stating a stability result for the continuous problem which was proved in [6].
(Ω) and X(t) ∈ H 1 (B) be solutions of Problem 1. Assume that ∂X(t)/∂t ∈ L 2 (B) and consider the elastic potential energy of the body given by Then the following conservation property is satisfied for almost every t ∈ (0, T )

Existence and uniqueness of the linearized problem
Not many results are available in the literature about existence and uniqueness of the solution to fluid-structure interaction problems. This is not surprising since the coupling between fluids and solids gives rise in general to highly non linear problems. In the case when a fluid is containing rigid solids or elastic bodies described by a finite number of modes, existence and uniqueness of weak solutions have been studied for instance in [22,25,26,27,30,36,37,39,47,48,49]; when a fluid is enclosed in a solid membrane then the existence and uniqueness of weak solutions have been discussed in [3,20,42,43]. Moreover, local-in-time existence and uniqueness of strong solutions for an elastic structure immersed in a fluid are proved in [23,24,46,16,17].
In this section we describe the analysis performed in [11] about the existence and the uniqueness of a linearization of Problem 1 in the case when Λ = H 1 (B) and the bilinear form c is equal to the scalar product in H 1 (B). This is a first step towards the analysis of the full problem which could make use of some fixed point strategy.
We consider a given function X that describes the motion of the solid. We assume that X belongs to C 1 ([0, T ]; W 1,∞ (B)), is invertible with Lipschitz inverse, and coincides with the identity at time t = 0, that is X(s, 0) = s. Moreover, we assume that the motion of the solid is compatible with the incompressibility constraint, that is det(∇ s X(t)) = 1 for all t.
We choose a linear model for the elasticity, namely P(F) = κF; moreover, we introduce a new variable w(t) equal to the velocity of the solid ∂X(t)/∂t, so that, after neglecting the convective term in the Navier-Stokes equation, we are led to the following problem.
Theorem 2. Under the assumptions reported above, there exists a unique solution to Problem 2 that satisfies the following regularity . The proof of this result is obtained by considering first a reduced problem where the unknowns p and λ are eliminated since the velocity is sought in the kernel of the divergence operator V 0 and the pair (u(t), w(t)) is required to satisfy the constraint In this setting, the proof follows a suitable modification of the Galerkin arguments used in [50] for the analysis of Navier-Stokes equations. Finally, the Lagrange multiplier and the pressure are recovered by using Lax-Milgram lemma and the Banach closed range theorem.

Time advancing schemes
We begin in this section the study of the numerical approximation of Problem 1, starting from the time discretization.
Let us introduce a time discretization parameter ∆t, and let us denote by t n , n = 0, . . . , N the corresponding nodes; the following system is obtained by the application of the backward Euler scheme.
where X −1 can be defined, for instance, from the following equation In [6] the following stability estimate was proved for the time discretization presented in Problem 3 Despite the nice stability property, it is clear that solving Problem 3 requires expensive numerical strategies in order to deal with the fully implicit non-linear scheme.
For that reason, we considered other semi-implicit schemes based on the use of the position of the structure at time n instead of n + 1. A possible semi-implicit version of (10) reads Moreover, in each particular situation, the quantity P(F n+1 ) should also need a linearization in order to avoid the presence of fully implicit terms.
The following stability estimate was proved in [6].
Proposition 3. Let us assume that the potential energy density W is a C 1 convex function, then the solution of (11) satisfies Remark 2. The stability results presented in Proposition 3 is a significant improvement over other schemes used for the approximation of fluid-structure interactions problems. A keystone result in this framework is reported in [19] where it is shown that schemes based on the Arbitrary Lagrangian Eulerian (ALE) approach cannot be stable, when the density of the fluid is close to that of the solid, unless they are fully implicit. Within the Immersed Boundary Method (IBM), when finite differences are used for the space discretization, it is shown that unconditional stability estimates can be obtained in some circumstances [44]. In our previous works we have shown a conditional stability, subject to a CFL condition, for the FE-IBM [7,12].
In [15] we investigated how to apply higher order schemes. We have to pay attention to the term involving the second time derivative of X; as it is common in this case, we reduce the order of the time derivative by introducing a new variable w corresponding to the first derivative of X (see also Problem 2). For instance, a scheme based on the BDF2 discretization reads (12) The following stability estimate was proved in [15].
Proposition 4. Let us assume that the Piola-Kirchhoff tensor is linear P(F) = κF, then the solution of (12) satisfies We refer the interested reader to [15] for other second order schemes based on the Crank-Nicolson method and to their corresponding stability properties which are analogous of the ones presented above. Some numerical experiments will be presented in Section 6.

Analysis and finite element approximation of the associated saddle point problem
In this section we discuss the finite element discretization in space of our problem. We consider the semi-implicit version of one of the schemes introduced in the previous section. At each time step we have to solve a stationary problem that we are going to present, approximate, and analyze. We consider X ∈ W 1,∞ (B) d that corresponds to X n and u ∈ L ∞ (Ω) that corresponds to u n . In this section we deal with the following Piola-Kirchhoff tensor Moreover, we define the following bilinear forms where the constants α, β, and γ depend on the time step and on the coefficients of our model. For instance, in the case of backward Euler method we have α = ρ f /∆t, β = δ ρ /∆t, γ = κ∆t.
Problem 4 (Saddle point problem). Given f ∈ L 2 (Ω), g ∈ L 2 (B), and d ∈ L 2 (B), find u ∈ H 1 0 (Ω), p ∈ L 2 0 (Ω), X ∈ H 1 (B), and λ ∈ Λ such that In general, f , g, and d are related to quantities at previous time steps. For instance, in the case of the backward Euler scheme we have Problem 4, after converting bilinear forms into linear operators with natural notation, reads  which has a saddle point structure. While the analysis of this problem has been published in [9], in [10] we observed that it was more convenient to rearrange the unknowns as follows  The saddle point structure is evident by introducing the following operators: A : V → V and B : V → L 2 0 (Ω) given by × Λ equipped with the graph norm. In particular the operator A has itself a saddle point structure which is highlighted by the dashed lines. In [10] it is shown that Problem 4 is well posed by proving the following properties.
• The operator A is invertible in the kernel of B.
• The operator B is surjective.
Since A is characterized by a saddle point structure, its invertibility is proved by showing the validity of two inf-sup conditions, while the surjectivity of B follows from the standard inf-sup condition of Stokes-like problems. For the sake of completeness, we recall the statements of the results that are needed in order to prove that A is invertible in the kernel of B.
We start by observing that V = (v, X, µ) ∈ V belongs to the kernel of B if and only if div v = 0. We recall that the divergence free subspace of H 1 0 (Ω) was denoted by V 0 .
In order to study the operator A we use the following kernel (see also (9)) and we show that there exists α 0 > 0 such that The invertibility of A in the kernel of B is then implied by the following inf-sup condition: there exists a constant β 0 > 0 such that The above estimate holds true for both choices of the bilinear form c defined in (6) and (7), and is a natural consequence of the definition of the norm of Λ.

Finite element discretization. The finite element discretization of Problem 4 is performed by considering finite dimensional subspaces
, and Λ h ⊂ Λ. We assume that the spaces V h and Q h are an inf-sup stable choice for the approximation of the Stokes problem.
In this paper we consider a more general setting than the one studied in [9] where we assumed that Λ h and S h were equal to each other.
The finite element spaces are constructed starting from three fixed shape-regular meshes: T V with mesh size h x for the domain Ω, T S with mesh size h s for the domain B, and T Λ with mesh size h λ for the domain B. The first mesh is associated with the use of the Eulerian variable x, while the other two meshes correspond to the Lagrangian variable s. Here we are assuming that Ω and B are polytopes and that B corresponds to the initial configuration of the solid. If this is not the case, then further approximations should be introduced. In any case a crucial property of our model is that the meshes are fixed during the entire evolution of the system.
The discrete counterpart of Problem 4 can be written as follows.
Problem 5 (Discrete saddle point problem). Given f ∈ L 2 (Ω), g ∈ L 2 (B), and In the formulation presented above we used the notationĉ for the discrete realization of the bilinear form c considered in Problem 4. Let us detail how this realization looks like in the two cases described in (6) and (7). If Λ = Λ 1 , observing that any reasonable finite element space Λ h is included in L 2 (B), it is possible to identify the duality pairing c 1 with the inner product of L 2 (B), so that we takê On the other hand, in the case Λ = Λ 2 we can take the same bilinear form as in the continuous caseĉ We are going to use the same notation c for both approaches as for the continuous case. When we need to refer explicitly to one of the two formulations, we shall use the full notation.
The analysis of the discrete problem makes use of the same technique that we described above for the continuous case. We report the main ingredients of the proof in a more general setting than it was presented in [10]; this is also the occasion to amend some detail of [9].
Using the notation introduced above for the space V, we introduce as follows the bilinear forms A : V × V → R and B : V × L 2 0 (Ω) → R in order to highlight the saddle point structure of the problem and to make easier the description of the result where we used the notation U = (u, X, λ) and V = (v, Y, µ). It is clear that the bilinear forms A and B correspond to the operators A and B defined above. We denote by V h = V h × S h × Λ h the subspace of V that we are using for the approximation. Hence, Problem 5 reads: given f ∈ L 2 (Ω), g ∈ L 2 (B), and Let V 0,h be the subset of V h containing the discretely divergence free vectorfields, that is v h ∈ V 0,h if and only if We are going to use the discrete kernel We state the following compatibility between the spaces S h and Λ h that will be useful in the sequel. Assumption 1. There exists a constant ζ > 0 such that for all µ h ∈ Λ h it holds (17) sup In order to show the stability of (16) we need to prove the following inf-sup conditions [4].
• There exists γ 1 > 0 such that (18) inf • There exists γ 2 > 0 such that The inf-sup condition for the bilinear form B is immediate if the spaces V h and Q h are a good Stokes pair. Indeed it is easy to see that where γ 2 is the inf-sup constant related to V h and Q h for the divergence operator.
In order to show the inf-sup condition for the bilinear form A , we start with the following proposition.
Proposition 5. For all β ≥ 0, there exists a constant α 1 > 0 not depending on the mesh sizes such that Proof. This proposition extends the conclusions of [10,Prop. 7]. For β > 0, the result follows directly from 0,B . The next step is to show that we can control X h 0,B by the right hand side of (20). This can be done at once for both possible choices of Λ and c. In order to use the Poincaré inequality we split X h as the sum of its mean valueX h and the rest, so that The constant partX h can be estimated by using the fact that the finite element space Λ h contains the global constant functions as follows.
Indeed, if µ h is constant then the term c(µ h , X h −X h ) vanishes and, even in the case when c is the scalar product in H 1 (B), the term involving ∇ s µ h vanishes so that c acts as the scalar product in L 2 (B). Hence, we get the final bound X h 0,B ≤ u h 0 . The next step consists in showing the following uniform inf-sup condition. Proposition 6. Let us suppose that Assumption 1 is satisfied. Then, for β 1 = ζ from (17) we have Proof. Using Assumption 1 we have We now present some possible choices of S h and Λ h for which Assumption 1 holds true.
We start by considering the case when Λ = Λ 2 and the bilinear formĉ = c 2 , namely it corresponds to the scalar product in H 1 (B). The most natural situation is when Λ h ⊆ S h which is the object of the following proposition. This condition is satisfied, for instance, if the mesh T S is the same as T Λ or a refinement of it and the space S h contains polynomials of degree higher than or equal to those in Λ h . Proposition 7. Let Λ = Λ 2 and the bilinear formĉ = c 2 be the scalar product in Hence the inf-sup condition (17) holds true with ζ = 1.
Let us now consider the case when Λ = Λ 1 is the dual of H 1 (B) and the bilinear formĉ =ĉ 1 is the scalar product in L 2 (B). We take again the most natural situation when Λ h ⊆ S h as in Proposition 7. In this case, however, the validity of the inf-sup condition (17) relies on an additional hypothesis that involves the H 1 (B)-stability of the L 2 (B)-projection onto S h . Proposition 8. Let Λ = Λ 1 = (H 1 (B)) and the bilinear formĉ =ĉ 1 be the scalar product in L 2 (B). Let P 0 denote the L 2 (B)-projection from H 1 (B) onto S h and assume that there is a constant C such that Then, if Λ h ⊆ S h , the inf-sup condition (17) is satisfied.
Proof. By definition of the norm in Λ there existsỸ ∈ H 1 (B) such that where in the last equality we used Λ h ⊆ S h . Finally, using the H 1 (B)-stability of P 0 stated in (21), we get Hence, the proposition is proved with ζ = 1/C 0 .
The cases considered in Propositions 7 and 8 generalize the situation discussed in [9], where Λ h was chosen equal to S h . It will be the object of further investigation to explore other possible combinations for Λ h and S h . In particular, it would be quite natural to take a space of discontinuous finite elements for the multiplier in the case when Λ = (H 1 (B)) . On the other hand, our present analysis does not cover for instance the situation when Λ h is the space of piecewise constants and S h is the space of continuous piecewise linear elements in each component: Assumption 1 requires dim(S h ) ≥ dim(Λ h ) as a necessary condition, which is not satisfied on general meshes for this choice of finite elements.
The results of this section can be summarized in the following stability and convergence theorems. If, moreover, V h and Q h satisfy the usual compatibility condition for the solution of the Stokes problem, then the inf-sup condition (19) holds true.
Proof. The results of this theorem follow from the previous propositions with classical arguments related to the stability of saddle point problems [4] (see also [51]).
The inf-sup condition (18) is the necessary and sufficient condition for the uniform invertibility of the matrix restricted to the discrete kernel of the matrix B = B f 0 0 , where the blocks A f , A s , B f , C f , and C s are matrix representations of the corresponding operators in (14).
Proposition 5 states the uniform invertibility of the block Proposition 6 states the surjectivity of this last matrix with uniform bound of its inverse.
Putting things together, we get the inf-sup condition (18). The second part of the theorem has been discussed after formula (19).
From the stability of the discrete problems, the convergence result follows in a straightforward way.

Numerical results
In this section we collect some numerical experiments that have been reported in previous papers and that confirm the effectiveness of the method.
We start with a test reported in [6] confirming the unconditional stability stated in Proposition 3. We consider a benchmark test problem where at the initial time the solid occupies an ellipsoidal region which evolves approaching a circular equilibrium configuration. We approximate the problem by using the enhanced Bercovier-Pironneau element introduced and analyzed in [5], consisting in a P1iso-P2 discretization of the velocities and in a continous P1 discretization of the pressures augmented by piecewise constant functions in order to improve the mass conservation of the scheme. We compare our fictitious domain approach FE-DLM (solid line) with the FE-IBM scheme (dashed line), see [12]. We take Ω equal to the square of side (−1, 1) and we study a ring-shaped immersed structure with reference configuration given by B = {x ∈ R 2 : 0.3 ≤ |x| ≤ 0.5}. For symmetry reasons, we reduce the computation to a quarter of Ω so that the configuration is the one reported schematically in Figure 1.
The materials properties are ρ f = 1, ν = 0.05, δ ρ = 0.3, and κ = 1. The solid mesh size is equal to 1/8 and the figures show the behavior of the following energy ratio as a function of the time step and of the fluid mesh size (22) Π(X n h , u n h ) = In Table 1 we report the results presented in [15] about the convergence rates in time when different time schemes are used. In these computations the mesh of Ω is based on a subdivision of (−1, 1) in 32 equal subintervals and the structure    (22)) for different ∆t when h x varies. The solid line corresponds to the formulation FE-DLM described in this paper, while the dashed line refers to the FE-IBM scheme which is only conditionally stable is modeled by a Lagrangian mesh obtained by halving the meshsize of the one reported in Figure 1. The fluid is initially at rest and the structure is stretched by a factor 1.4 in the vertical direction and shrunk by the same factor in the horizontal direction.
The physical parameters are ρ f = ρ s = 1, ν = 0.1, κ = 10, and T = 1. We consider BDF1 (semi-implicit backward Euler (11)), BDF2 (see (12)), and two variants of Crank-Nicolson scheme. We denote by CNm the case when the nonlinear terms are evaluated using the midpoint rule and by CNt the case when the trapezoidal rule is used.
The reference solution is calculated by using a smaller timestep with the BDF2 scheme.