Time-Continuous and Time-Discontinuous Space-Time Finite Elements for Advection-Diffusion Problems

We construct four variants of space-time finite element discretizations based on linear tensor-product and simplex-type finite elements. The resulting discretizations are continuous in space, and continuous or discontinuous in time. In a first test run, all four methods are applied to a linear scalar advection-diffusion model problem. Then, the convergence properties of the time-discontinuous space-time finite element discretizations are studied in numerical experiments. Advection velocity and diffusion coefficient are varied, such that the parabolic case of pure diffusion (heat equation), as well as, the hyperbolic case of pure advection (transport equation) are included in the study. For each model parameter set, the L2 error at the final time is computed for spatial and temporal element lengths ranging over several orders of magnitude to allow for an individual evaluation of the methods' spatial, temporal, and spacetime accuracy. In the parabolic case, particular attention is paid to the influence of time-dependent boundary conditions. Key findings include a spatial accuracy of second order and a temporal accuracy between second and third order. The temporal accuracy tends towards third order depending on how advection-dominated the test case is, on the choice of the specific discretization method, and on the time-(in)dependence and treatment of the boundary conditions. Additionally, the potential of time-continuous simplex space-time finite elements for heat flux computations is demonstrated with a piston ring pack test case.


Motivation
Multiple features make space-time finite elements an attractive solution strategy for time-dependent partial differential equations (PDE).First, space-time finite elements provide a uniform framework for error analysis as no distinction is made between spatial and temporal coordinates 1 , which can also be used in adaptive refinement of the combined space-time mesh 2 .Moreover, spacetime finite elements allow for parallel-in-time (PinT) computations which have inherently more potential for parallelization than spatial finite elements combined with a sequential time-stepping scheme 3 .Furthermore, space-time finite elements are a natural 0 Preprint submitted for publication arXiv:2206.01423v1[math.NA] 3 Jun 2022 choice to discretize time-dependent spatial computational domains, e.g., in fluid-structure interaction (FSI) simulations 4,5,6,7 .In particular, simplex space-time finite elements 8 can provide a boundary conforming space-time mesh for spatial domains that change topology over time 9 .To benefit from these advantages, space-time finite elements have been used to perform simulations in various fields of computational fluid dynamics (CFD).Recent examples of simplex space-time simulations include the computation of complex fluid flows in production engineering applications 10,11 and the computation of dense granular flows 12 .Likewise, compressible flows have been successfully simulated on unstructured space-time meshes 13,14,15 .Note that the solution of transient three-dimensional problems with space-time finite elements requires four-dimensional meshes.Recent advances in generation 9,16,17 , adaptation 18 , and numerical handling 19,20 of four-dimensional meshes mark the state-of-the-art in this active research field.For efficiency considerations and refinement strategies, it is important to know the convergence behavior of the space-time finite element solution towards the physical or analytical solution of the simulated test case.However, for simulations based on the incompressible or compressible Naiver-Stokes equations it is an intricate task to estimate exact convergence orders, since numerical reference solutions can be influenced by round-off errors or implementation issues.Instead, we consider in this paper advection-diffusion problems -which lend themselves to an analytical solution -as a prototype for more complex flow problems 21 .Based on the advection-diffusion equation, one can investigate the performance of numerical schemes with respect to transient, advective, and diffusive effects as well as their interplay.Besides, advection-diffusion equations also model a variety of physical problems, e.g., the concentration of a chemical species transported by an ambient flow or the temperature of a fluid streaming along a heated wall 21 .Therefore, it is of great interest to analyze the convergence behavior of numerical schemes for advection-diffusion problems.

Literature Review
Shakib and Hughes 22 present a Fourier analysis of space-finite elements with tensor-product structure applied to an advectivediffusive model problem with periodic boundary conditions.The method is found to be third order accurate with respect to the time step size for the pure advection and pure diffusion case.A summary of space-time finite element methods for convective transport problems is provided by Donea and Huerta along with numerical tests 23 .Moreover, linear tensor-product space-time finite elements can be related to a spatial discretization with finite elements and a temporal discretization with the Crank-Nicolson scheme 8 .Studies of this resulting method often focus either on parabolic problems (heat equation) 24 or on the pure advection case (transport equation) 25 .Moreover, a Crank-Nicolson type space-time finite element method for evolution problems on moving meshes is proposed and analyzed by Hansbo 26 .The method uses tensor product elements that are inclined in space-time with a slope given by the convection velocity.It is reported that the aligned space-time orientation improves the precision and facilitates the solution of the discrete system.Focusing on the parabolic limit case, time-continuous tensor-product space-time finite elements have been analyzed by Aziz and Monk 27 .In more recent works, also unstructured space-time finite elements which do not require any tensor-product structure are addressed, e.g., by Steinbach 28 .Furthermore, Langer and Schafelner2,29 investigate the scaling behavior of unstructured spacetime finite element methods for parabolic problems in parallel computations.Note that this work is also extended to hexahedral space-time discretizations 30 .Moreover, Langer and Zank propose and investigate new efficient direct solvers for time-continuous tensor-product discretizations of the parabolic initial boundary value problem 31 .The influence of linear constraints, e.g., timedependent Dirichlet boundary conditions, on discontinuous Galerkin time discretization methods for parabolic problems is treated by Voulis and Reusken 32 .

Scientific Novelty and Limitations
To the best of the authors' knowledge, there is no previous comprehensive numerical study that analyses the convergence behavior of tensor-product and simplex-type finite elements for the complete range of model parameters of advection-diffusion problems and for spatial and temporal element sizes over several orders of magnitude.On the one hand, the computational evaluation of the convergence behavior is advantageous in the sense that a simple variation of the input parameters allows to switch from a parabolic to a hyperbolic problem.Therefore, the computational approach facilitates a study of the precise influence of parameter variations.On the other hand, the numerical study is limited to specific test cases and for those considers only the

Paper Organization
In the remainder, we proceed as follows.In Section 2, four space-time discretizations are presented and descriptive naming is proposed.In Section 3, we apply the methods to an initial boundary value problem based on the advection-diffusion equation.Section 4 collects the results of a computational error analysis of the time-discontinuous discretizations and compares the results with the theoretically expected convergence behavior.In Section 5, we demonstrate the particular potential of simulations on time-continuous simplex space-time meshes in a piston ring pack application.Concluding remarks are offered in Section 6.

METHOD CLASSIFICATION
To introduce the specific space-time discretizations investigated in this work, the naming of involved entities is briefly reviewed 1,15 .We consider a spatial computational domain Ω ⊂ ℝ sd , where sd denotes the number of spatial dimensions.That domain Ω and a time interval, = [0, ] ⊂ ℝ, span the space-time continuum ⊂ ℝ sd +1 .In the following, we consider four ways to approximate the solution of PDEs on with finite elements.Sample slicings ℎ of the space-time domain = [ 0 , 3 ] × [ 0 , 3 ] are shown in Figure 1.For the sake of clarity, the spatial domain Ω remains constant over time in these drawings.However, the proposed methods can also be applied to time-dependent spatial domains Ω( ) 4,9,16 .The first two discretization techniques (Figure 1a and 1b) seek an approximation that is continuous across .In contrast, the second two (Figure 1c and 1d) seek an approximation that is discontinuous at certain times, which leads to a discontinuous Galerkin method for the temporal discretization.In these time-discontinuous cases, is sliced into space-time slabs .As indicated in the drawings of Figure 1c and 1d, the boundary of each space-time discretization consists of three parts: the spatial discretization at the lower time level Ω ℎ l = Ω ℎ ( = ), the spatial discretization at the upper time level Ω ℎ u = Ω ℎ ( = +1 ), and the discretization of   the space-time boundary ⊂ ℝ sd which is the temporal evolution of the spatial domain boundary Γ ⊂ ℝ sd −1 .The size of spacetime slabs in temporal direction is denoted by Δ .To later apply one uniform finite element formulation for the time-continuous and time-discontinuous cases, we regard the complete space-time domain in the time-continuous case as space-time slab 0 .Both time-discretization approaches can be combined with prismatic elements with tensor-product structure, or simplex elements.The combinations form the four discretization methods C-PST, C-SST, D-PST, and D-SST.In PST methods, a discretization of with prismatic space-time elements can be easily obtained by extrusion of a spatial discretization of Ω in time.C-PST is a continuous finite element discretization in space and time as described by Aziz and Monk 27 .When combined with linear shape functions, it is also known as cg(1)cg (1).The time-discontinuous D-PST method is also referred to as cg(1)dg(1) for example by Quarteroni et al. 33 .SST discretizations can be generated by subdividing the prismatic elements into simplex elements (Figure 1b and 1d).More complex SST mesh generation procedures also allow for local temporal refinement by node insertion 8 or fully-unstructured space-time meshes 15 as shown in Figure 2.For each space-time slab , an 1 -conformal finite element approximation space 1 ℎ, is constructed based on one of the discussed discretizations and element basis functions 15 .In case of PST discretizations, we consider the ℙℝ 1 and ℚ 1 basis functions of the simplex-based prismatic and cuboid Lagrange finite element.In case of SST discretizations, we use the ℙ 1 basis functions of the simplical Lagrange finite element 34 .The four space-time discretizations introduced, are now employed in the solution of advection-diffusion problems.

APPLICATION TO ADVECTION-DIFFUSION EQUATION
We consider the time-dependent linear advection-diffusion equation Therein, the scalar unknown, ( , ) is a function of the spatial coordinates ( = ( , , ) for sd = 3) and time.The advection velocity is a given vector , and the diffusion coefficient is denoted by .As usual, the Laplacian of abbreviates Δ = ∇ ⋅ ∇ , based on the spatial gradient ∇ .Advection velocity and diffusion coefficient can be varied, such that the parabolic case of pure diffusion ( = ), as well as the hyperbolic case of pure advection ( = 0) are included.In the former case, Equation (1) is the heat equation, in the latter case the transport equation.Furthermore, the above equation lends itself to an analytic solution, hence, facilitating a computational error analysis as presented in Section 4.
A general characterization of advection-diffusion problems can be achieved with the dimensionless Péclet number Therein, a scalar measure of the advection speed, = ‖ ‖, is related to the diffusion coefficient and scaled by a characteristic length .As the dimensionless number compares the importance of advective and diffusive effects for a given test case, one can typically expect solutions with smaller gradients for test cases with lower Péclet number (when diffusion dominates).
To construct an initial boundary value problem, let us consider again a computational domain as for example shown in Figure 2b.The associated space-time boundary is assumed to consist of a Dirichlet part and a Neumann part , such that = ∪ and ∩ = ∅.Then, we obtain an initial boundary value problem, as we require Equation (1) to hold on , along with a known initial condition 0 and given Dirichlet boundary conditions .The complete statement of the initial boundary value problem reads ( When applying one of the discretization techniques described in Section 2 to , the initial condition is enforced on Ω of the space-time slab 0 .The part of a space-time slab , where Dirichlet boundary conditions are prescribed is denoted by .A suitable interpolation of the Dirichlet boundary data ℎ allows us to define the trial function space and the test function space Considering that for time-discontinuous discretization methods the finite element approximation is discontinuous at the spacetime slab boundaries Ω and Ω , let ℎ ± abbreviate lim →0 ℎ ( ± ).Using these definitions, a discretized weak form of the initial boundary value problem can be stated as follows: For given initial conditions ℎ − 0 = ℎ 0 , find ℎ ∈  ℎ, such that on each time slab and for all ℎ ∈  ℎ, In the weak form above, the diffusion term was modified using integration by parts.The resulting boundary integral vanishes, since the test functions vanish on and homogeneous Neumann boundary conditions are assumed on .Moreover, the initial condition as well as the continuity of ℎ between time slabs is weakly enforced with the integral over the spatial computational domain Ω , the so-called jump term.The stability of the formulation is achieved with a SUPG term in the fourth integral 35 .We define the stabilization parameter SUPG as which accounts for local characteristics of the initial boundary value problem.In the first term, the space-time element metric is used to include directional element length information.The metric tensor, is based on the inverse of the Jacobian associated with the mapping from reference coordinates, ( , ), to physical coordinates, ( , ).Moreover, the metric tensor includes a square matrix of size sd +1, which accounts for the mapping to a regular reference element counteracting the influence of the element's node numbering 15 .A further analysis of node-numbering invariant element length measures for simplex elements is presented by Takizawa et al. 19 .Explicit forms of for simplex elements read for For other element types, an appropriate matrix is substituted.Recalling that accounts for the mapping to a regular reference element, it is clear that discretizations with pure tensor-product reference elements (ℚ 1 ) do not need an additional mapping-as the reference element is already regular.Therefore, can simply be replaced by the identity matrix.In case of a simplex-based prismatic reference element (ℙℝ 1 ), the partial tensor-product structure of the reference element is reflected in the choice of as shown below In the second term of Equation ( 7), the diffusive contribution to SUPG requires a measure of the spatial element length ℎ .For all considered element types, the length ℎ is obtained from the spatial part of the metric tensor = [ ] sd × sd as where the colon operator denotes the double contraction ∶ = ∑ , ⋅ .Moreover, the constant inv scales the diffusive contribution to SUPG .Inspired by an inverse estimate inequality proven in 36 , we chose for ℙ 1 and ℙℝ 1 discretizations For ℚ 1 discretizations, we use inv ≈ 1.To improve the consistency of our formulation in combination with linear finite elements, the second-order derivatives in the residual res( ℎ ) are obtained with a least-squares recovery technique 37 .For the parabolic case ( = ) and linear approximation functions, the weak form in Equation ( 6) is very close to the locally stabilized space-time finite element method presented by Langer and Schafelner in 29 Section 3.Only the definition of the stabilization parameter, SUPG or Θ , 2 Remark 13.4, and the enforcement of the initial condition differ.To provide a first test case for the four space-time discretization methods, we analyze the transient one-dimensional model problem We consider a time interval = ]0, 2] and the spatial computational domain Ω spans from −1 to 1.The model problem is characterized by the periodic boundary conditions and has the analytical solution The test case setup of IBVP 1 is also discussed by Mojtabi and Deville 38 and on a shifted computational domain by Shakib and Hughes 22 .
In the numerical solution procedure, we discretize the computational domain with eight elements in spatial and temporal direction as shown in Figure 3. Due to the periodic boundary conditions (−1, ) = (1, ), this leads to eight independent degrees of freedom in spatial direction.As the initial condition is enforced weakly, the time-continuous discretizations have nine nodes in time direction with one degree of freedom each.Therefore, C-PST and C-SST simulations use 9 × 8 = 72 degrees of freedom   Comparing the solution of D-PST in Figure 3c with the C-PST solution in Figure 3a, one can note jumps in the solution at the interfaces between the space-time slabs.These small discontinuities in the solution are in line with the weak enforcement of the continuity requirement in the weak form (Equation ( 6)).Also the D-SST solution is discontinuous at the interfaces between space-time slabs.However, these jumps are less pronounced and not visible in the rendering of Figure 3d.Regarding the SST discretizations (in Figure 3b and 3d), we can note that the solution is advected along the diagonal edges of the SST discretizations.
In this particular case with Δ = Δ and = 1, the characteristics perfectly align with the finite element edges.Figure 4 compares the numerical solutions ℎ of the four space-time discretization methods to the analytical solution at the final time .In the plot of the differences − ℎ (Figure 4a), the interpolation error between the nodal values is very prominent.Please, note that this error is inherent to the linear interpolation of a trigonometric function.Removing this unavoidable error (for linear approximation functions), Figure 4b connects the values at the finite element nodes with straight line segments.For the employed, very coarse discretizations, the nodal differences of the SST solutions to the analytical solution are smaller, despite the smaller number of degrees of freedom in comparison to the PST methods.Additionally, the D-SST method shows hardly any phase error.Comparing Figure 4a and Figure 4b, one can observe that the error of the finite element solution at the nodes is of the same order as the interpolation error.
Returning to the complete space-time solution (Figure 3), all four space-time discretization methods arrive at similar results.Given the extremely coarse discretization, one can consider all numerical solutions to be in accordance with the analytical solution.We therefore conclude that all four space-time discretization schemes are suitable for advection-diffusion problems.
The C-PST method has been analyzed for the heat equation theoretically and with numerical experiments by Aziz and Monk 27 .It is found that the use of linear finite element approximation functions in C-PST leads to a version of the Crank-Nicolson method.Moreover, the tensor-product approach of C-PST leads to a global linear equation system with specific structure for parabolic initial boundary value problems.This structure can be exploited in the construction of an efficient parallel solver as shown by Langer and Zank 31 .Still, we will not further consider the scheme in this paper.The C-SST method allows for space-time adaptivity on unstructured meshes 29 and in Section 5 the C-SST method is used to include topology changes of the spatial computational domain Ω in a boundary-conforming space-time mesh.However, Section 4 focuses on the time-discontinuous methods, D-PST and D-SST.

COMPUTATIONAL ERROR ANALYSIS OF TIME-DISCONTINUOUS DISCRETIZATIONS
To investigate the convergence of the time-discontinuous space-time discretizations D-PST and D-SST, a computational error analysis is performed.The following space-time convergence studies consider two test cases.Before investigating IBVP 1 (Equation ( 13)) for six model parameter sets in Section 4.2, we first consider the parabolic case ( = 0) of a second initial boundary value problem IBVP 2 with time-dependent Dirichlet boundary conditions in Section 4.1.Both initial boundary value problems have analytical solutions, which serve particularly well as reference solutions in the convergence studies, since they are independent of implementation issues or round-off errors introduced in computer arithmetic.
For each model problem, parameter set, and discretization method (D-PST, D-SST) a space-time convergence study with 198 simulations is performed.The numerical simulation settings are obtained as follows.We divide the computational domain in time direction (up to the final time = 2) into ts space-time slabs of constant size Δ = ∕ ts .We consider 15 levels of recursive temporal refinement such that ts is doubled from the coarser to the finer level In the same manner, the spatial domain is divided into ex elements of constant size Δ = 2∕ ex .The number of elements in spatial direction is given by For each simulation, the relative 2 error at the final time = 2 for the spatial refinement level and temporal refinement level is evaluated.In practice, we use an element-wise two-point Gaussian quadrature for the spatial integration Additionally, we measure the nodal errors as with the index running over all nodes, except for the last one.In model problem IBVP 1, node 1 and node ex + 1 have identical solution values enforced by the periodic boundary conditions.For both model problems, the prefactor with the 2 norm of the solution at the final time = 2 reads Simulations are performed for whole-numbered parameter pairs ( , ) corresponding to grid line intersections in Figure 5.To avoid unnecessary computational cost, we omit combinations of the finest refinement levels as shown in Figure 5.Note that the patch color is based on the mean value of the 2 error of the four simulations connected by a patch.To check for spatial convergence, we consider the finest temporal refinement level = 18 and vary = 4, … , 12, which corresponds to the bottom line of the plot in Figure 5. Analogously to investigate temporal convergence, we consider the finest spatial refinement level = 18 and vary = 4, … , 12.This corresponds to the rightmost line of the plot in Figure 5. On the space-time diagonal = , the numerical values of Δ and Δ coincide.Despite the different units that one would assign to the physical quantities, we use Δ = Δ to express that the numerical values are equal.Along the curve Δ = Δ , twelve data points = = 4, … , 15 are generated.

Parabolic model problem IBVP 2
In this section, we study the pure diffusion case of the model problem ( , ) = cos( ), at = 0, With the time-dependent Dirichlet boundary conditions ( ), IBVP 2 has the analytical solution The considered advection-diffusion equation ( 1) reduces for = 0 to the heat equation.For a corresponding initial boundary value problem with homogeneous boundary conditions, convergence estimates for D-PST discretizations are known from literature.Thomée presents in 39 Theorem 12.7 a superconvergence result for the temporal discretization error at the final time .
Considering linear basis functions, the error bound for the parabolic problem can be summarized as where is a positive constant independent of Δ and Δ .
In the following, we compare our computational findings to the theoretical result above.The results of the space-time convergence studies are visualized in convergence surfaces (see Figure 6).The surfaces are obtained by plotting the 2 error in logarithmic scale over the spatial and temporal refinement level indices and .Corresponding convergence surfaces based on the nodal error measure can be found in the Appendix A in Figure A1.For both discretization methods, the error plots result in a continuous surface (Figure 6a and 6b).Moreover, the surfaces show, that the error values in the area of the diagonal ( = , Δ = Δ ) are influenced by the spatial and temporal mesh size.However, on the finest spatial discretization level ( = 18), the error varies only with Δ .The same holds for the finest temporal refinement level ( = 18) and Δ .Therefore, extracting the curves = 18 or = 18 from the convergence surfaces gives us the isolated spatial or temporal convergence behavior of the methods.At first, focusing on the spatial convergence rates shown in Figure 7a, we observe a second-order spatial accuracy for both methods as the curves of D-PST and D-SST coincide.This is to be expected, as the same spatial mesh is used.Moreover, this observation is also in line with the theoretical result given in Equation (22).Next, looking at the temporal convergence rates in Figure 7b, a second order temporal convergence is observed for both methods.This is in strict contrast to the third-order time accuracy expected from Equation (22).
As pointed out by Voulis and Reusken in 32 , the reduced convergence order is due to the time-dependent boundary conditions.Moreover, it is shown in their work that superconvergence can be recovered by applying a temporal interpolation operator to the boundary condition ( ).The use of this interpolated boundary condition is equivalent to the time-discontinuous discretization of the boundary condition ( , ) = ( ).In our considered test case, the temporal convergence can be improved with the following treatment.On the upper time level Ω = Ω( ) of each space-time slab, the boundary condition is precisely evaluated as On the lower time level Ω , a modified boundary condition ̃ ( ) is applied.The modified boundary condition is constructed such that the linear interpolation of the finite element shape functions leads to the correct analytical mean of the boundary condition For the considered example, this yields Repeating the space-time convergence study with modified boundary conditions, we obtain the results shown in Figure 8.Here, D-PST reaches significantly smaller error values in comparison to the case shown in Figure 6a.Returning in the line plots of Figure 7a, it can be seen that the second order spatial convergence of both methods is not affected by the boundary condition treatment as all four curves coincide.But, for D-PST with ̃ ( ), third-order temporal convergence is indeed obtained (Figure 7b).The result numerically confirms that temporal superconvergence (as stated in Equation ( 22)) can also be obtained for timedependent boundary conditions with a proper treatment 32 .For D-SST with ̃ ( ), only quadratic temporal convergence is observed.The lower convergence order of the D-SST method with treatment of the time-dependent boundary conditions hints at the fact that superconvergence of the D-PST method is linked to the tensor-product structure of the discretization.However, also in case of the D-SST discretization, the proposed treatment of time-dependent boundary conditions is helpful-the error values decrease by approximately 25%.The purely spatial or temporal refinements are interesting as they show an isolated spatial or temporal convergence behavior, but they are certainly not efficient in terms of computational cost that is required to obtain a certain level of accuracy.Elaborating on this, we assume that the computational cost of a simulation is related to the number of degrees of freedom dof .For the considered discretizations, dof can be expressed by the number of time steps ts and the number or elements in -direction ex as Since we estimate the total computational cost by the number of degrees of freedom, it is of the order ( ).This can be used to balance the spatial and temporal discretization to minimize the computational cost for a desired error.The optimal relation between the spatial mesh size and the temporal therefore depends on the relation between the spatial and the temporal convergence order.If the spatial and temporal convergence order match, then the choice Δ = Δ is optimal.However, if we consider the setting in Equation ( 22), then the optimal choice is Δ 2 = Δ 3 .In the visualizations of the convergence studies, e.g., in Figure 8, we can identify the best space-time refinement strategy as the steepest decent in the convergence surfaces.In Figure 8b, an advantageous space-time refinement strategy for D-SST essentially follows Δ = 4 ⋅ Δ .In contrast, for D-PST with ̃ ( ) (Figure 8a), the second-order spatial accuracy and third-order temporal accuracy lead to an advantageous space-time refinement strategy along the curve Δ 3 = Δ 2 .The 2 errors along the curves Δ = 4 ⋅ Δ and Δ 3 = Δ 2 are extracted from the convergence surfaces and plotted in Figure 9. On the space-time diagonal with offset (Δ = 4 ⋅ Δ , Figure 9a), both methods show a second-order convergence for the computations with and without ̃ ( ).For D-PST with ̃ ( ), the curve lies is the zone where the spatial error dominates, hence, we expect second-order convergence also for this method.However, in contrast to Figure 7a, the curves do not coincide and the treatment of the time-dependent boundary conditions proves advantageous in terms of the absolute error values.Note that transitioning from one data point to the next along the space-time diagonal doubles ts and ex .Following the advantageous refinement strategy for D-PST with ̃ ( ), Figure 9b shows that the third-order temporal accuracy of D-PST with ̃ ( ) is retained along the curve Δ 3 = Δ 2 .Summarizing the parabolic model problem analysis, both methods, D-PST and D-SST, converge at least quadratically against the analytical solution.With proper treatment, D-PST converges cubically with respect to Δ even for time-dependent boundary conditions.

Advective-diffusive model problem IBVP 1
While the main challenge in the previous Section 4.1 was the treatment of time-dependent boundary conditions, this section investigates the convergence behavior of the methods as the model parameters transition from the parabolic case to advectiondiffusion cases and to the hyperbolic case.The numerical error analysis of D-PST and D-SST is therefore continued with the model problem IBVP 1 (Equation 13).Six parameter sets are considered.They include the parabolic case = 0, = 0.1, Pe = 0,
next to four advection-diffusion cases with decreasing viscosity = 1, = 0.1, 0.01, 0.001, 0.0001, Pe = 1, 10, 100, 1000 and the hyperbolic case = 1, = 0, Pe = ∞.The periodic boundary conditions do not require the treatment of time-dependent boundary conditions.As before, the results of the space-time convergence studies are visualized in convergence surfaces.Four representative convergence surfaces are shown in Figure 10.The complete set of twelve surfaces can be found in Figure A2 in the Appendix.For all parameter sets, a continuous surface is obtained.Furthermore, for both space-time discretizations, the advection-diffusion cases with increasing Péclet number present a smooth transition from the pure diffusion to the pure advection case.Comparing the parabolic cases (Figure 10a and Figure 10b), the convergence surfaces of D-PST and D-SST clearly differ.D-PST reaches smaller error values due to the superconvergence of the discretization with tensor-product elements (Equation ( 22)).The hyperbolic cases (Figure 10c and Figure 10d) show only a slight difference for the simulations on coarse meshes with Δ = Δ , yet, on the finer meshes D-PST and D-SST arrive at very similar results.After this brief view on the convergence surfaces, we now analyze spatial and temporal convergence by means of line plots.Spatial convergence results are presented in Figure 11.Both space-time methods converge quadratically with respect to Δ for the complete model parameter range from Pe = 0 up to Pe = ∞ and for all values of Δ .Besides the constant convergence rates, there is an influence of the Péclet number on the actual relative error values.The solutions for the more diffusive cases, are slightly more accurate.Figure 12 collects the results of the temporal convergence study.For this model problem, D-PST is observed to converge cubically with respect to the time step for the complete range of Péclet numbers, despite the fact that the curve of Pe = 0 is shifted to smaller error values as shown in Figure 12a.This behavior is in line with the results obtained by Shakib and Hughes in a Fourier analysis of the purely advective and purely diffusive limiting case of this model problem 22 .Note that the specific mesh connectivity (stencil) of the D-PST discretization is used in the Fourier analysis and the results hence do not apply to a D-SST discretization.The D-SST results presented in Figure 12b show a strong influence of the Péclet number.The method is second-order time accurate in the parabolic case and third-order accurate in the hyperbolic case.For advection-diffusion cases, we observe a smooth transition of the convergence behavior from second to third order.However, rather than converging at a constant intermediate rate, D-SST converges for the advection-diffusion cases cubically up to some Δ turn (Pe), where the convergence rate transitions to two.For smaller Péclet numbers, the transition occurs at larger time steps, which is earlier in the convergence history.So far, we discussed the convergence of the 2 errors, but, also the convergence surfaces of the nodal error measure (Equation ( 18)) show interesting features of the discretization methods.The complete set of the twelve nodal error based convergence surfaces can be found in Figure A3 in the Appendix A. The results of the hyperbolic case are presented in Figure 13.The nodal error visualization of the D-PST results (Figure 13a) is a continuous surface as for the 2 error.In contrast, the D-SST results (Figure 13b) show a strong discontinuity for the simulations with Δ = Δ .For these cases, the space-time finite element edges align with the characteristic curves along which the solution is transported.We observe that the finite element approximation coincides with the exact solution at the nodes (up to a round off error < 1.0 × 10 −10 ) for all refinement levels.This astonishing behavior is described by Demkowicz and Oden as 'extra superconvergence' 40 .Away from the diagonal = , error values are obtained that are similar to the ones of D-PST.Spatial convergence results of the D-PST method in nodal error measure are extracted as line plots and shown in Figure 14a.We see once more a strong influence of the Péclet number on the convergence behavior.In the pure advection case, the method converges in the nodal error measure with fourth order up to a Δ turn (Pe) and then transitions to second order.For smaller Péclet numbers, the transition occurs at larger Δ , i. e., earlier in the convergence history.In the simulations with the fourth order convergence relative to each other, the element size in time direction Δ is very small compared to Δ .In consequence, the small Δ leads to such a small stabilization parameter SUPG (Equation ( 7)), that the influence of the SUPG term vanishes and the Galerkin method is recovered.The nodal error of the Galerkin method for the pure advection case is fourth order accurate with respect to Δ , as shown in the Fourier analysis of Shakib and Hughes 22 .In the pure diffusion case, the method is second order accurate with respect to Δ over the entire element size range.Figure 14b shows nodal error results of the D-SST method for the six model parameter sets along the space-time diagonal Δ = Δ .Most notable is the 'extra superconvergence' of the SST method with characteristics aligned element edges for the pure advection case.For the five other model parameter sets, we observe a second order space-time convergence along the diagonal Unfortunately, it is highly unlikely that the finite element edges of higher-dimensional space-time meshes are aligned with the solution characteristics for general flow conditions.Therefore, we come to the following outlook for problems of engineering interest.Under the assumption that our findings carry over from the scalar one-dimensional advection-diffusion cases to higherdimensional cases modeled with (in)compressible Navier-Stokes equations, we expect a spatial accuracy of second order and a temporal accuracy between second and third order.As we have seen, the temporal accuracy of the time-discontinuous spacetime methods tends towards third-order depending on how advection-dominated the test case is, on the element type used for discretization, and on the time-(in)dependence and treatment of the boundary conditions.(ii), and to dissipate heat from the piston to the surrounding cylinder to prevent overheating of the piston (iii).In this test case, we investigate the heat flux in a simplified model of a piston ring pack .
Figure 16 shows the two-dimensional geometry of a schematic piston ring pack with only one ring.The considered geometry includes a part of the piston around the groove in which the piston ring is located, as well as a part of the cylinder liner which comes into contact with the piston ring.As shown in Figure 16, the piston ring is represented by a square with a generic side length of 0.5; its corners are rounded with radius = 0.1.In the following simulation, we investigate the conductive heat transfer in the metal parts and across the contact interfaces between piston, ring, and liner.The heat transfer is modeled with the parabolic case of Equation (1), i.e., the advection velocity is set to zero and we obtain a Péclet number of zero.The thermal diffusivity in the solids = , that accounts for the thermal conductivity , density , and specific heat , is here modeled with a generic diffusion coefficient = 0.495, as employed in Equation ( 1).The test case is further characterized by the temperatures where ℎ( ) denotes the Heaviside function.
What makes this test case challenging is the ring motion.In the course of an engine working cycle, the piston ring is in contact with different parts of the piston and the liner.We consider a prescribed ring motion defined by the ring center position, ( ( ), ( )), as shown in Figure 17.During the simulated time interval ∈ (0, 2.6), the ring is first in contact with the upper edge of the piston groove, then moves downwards and is free-floating for ∈ (0.2, 0.6), before it comes into contact with the lower edge of the piston groove.These three states are also visualized in the first three figures in the left column of Figure 18.
In the following, the ring moves towards the liner, slides upwards along the liner and finally returns to the initial position.In a C-SST approach 9 , the given ring motion is included in the computational space-time domain as shown in Figure 15.We used GMSH 41 to discretize the domain with a fully unstructured space-time mesh.The resulting mesh consists of 151,911 tetrahedral elements connecting 35,341 nodes.Moreover, the mesh is refined in areas where large spatial and temporal solution gradient are expected, i.e., the curves where the ring comes into contact with the piston and the liner in the course of the simulation.The simulation results are collected in Figure 18 and Figure 19. Figure 18 shows the temperature distribution in the piston ring pack at eight time instances.Most of the time, the temperature solution in the piston and liner parts closely follows the prescribed boundary conditions.Larger spatial temperature variations are primarily encountered in ring.In particular at = 2.0, the ring directly connects the hot upper groove edge of the piston with the cooler liner.As indicated by the large temperature gradients, this configuration leads to the maximal conductive heat transfer.In Figure 19, the temperature at the piston ring center, ( ( ), ( )), is plotted over time.It is observed that the temperature is approximately constant during the interval ∈ (0.2, 0.6), which is expected as there is no conductive heat transfer to or from the free-floating ring.The strongest decrease in temperature is observed during the interval ∈ (1.2, 1.4), when the ring first comes Time-continuous prismatic space-time method (C-PST).
Local temporal refinement in D-SST mesh.

FIGURE 2
FIGURE 2Variants of simplex space-time discretization methods.

FIGURE 3
FIGURE 3 Solution of IBVP 1 computed with four space-time discretization methods.

( a )FIGURE 4
FIGURE 4 Comparison of four space-time discretization methods for IBVP 1.

1 FIGURE 5
FIGURE 5 Organization of space-time convergence study based on 2 error .

FIGURE 6
FIGURE 6 Convergence visualization of 2 error for parabolic problem configuration IBVP 2.

FIGURE 8
FIGURE 8 Convergence visualization of 2 error for IBVP 2 with modified boundary condition.

FIGURE 10
FIGURE 10 Convergence visualization of 2 error for model problem IBVP 1.

FIGURE 12
FIGURE 12Temporal convergence for IBVP 1 for six model parameter sets.

6 FIGURE 18 FIGURE 19
FIGURE18 Temporal evolution of temperature in piston ring pack .

FIGURE A1
FIGURE A1Convergence of nodal error measure for parabolic problem IBVP 2.

FIGURE A2
FIGURE A2Convergence visualization of 2 error for model problem IBVP 1.

FIGURE A3
FIGURE A3Convergence visualization of nodal error for model problem IBVP 1.